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Abstract 

This is a written version of a series of lectures aimed at graduate students in 
astrophysics/particle theory/particle experiment. In the first part, we explain 
the important progress made in recent years towards understanding the ex- 
perimental data on cosmic rays with energies > 10 8 GeV. We begin with a 
brief survey of the available data, including a description of the energy spec- 
trum, mass composition, and arrival directions. At this point we also give a 
short overview of experimental techniques. After that, we introduce the fun- 
damentals of acceleration and propagation in order to discuss the conjectured 
nearby cosmic ray sources, and emphasize some of the prospects for a new 
(multi-particle) astronomy. Next, we survey the state of the art regarding the 
ultrahigh energy cosmic neutrinos which should be produced in association 
with the observed cosmic rays. In the second part, we summarize the phe- 
nomenology of cosmic ray air showers. We explain the hadronic interaction 
models used to extrapolate results from collider data to ultrahigh energies, and 
describe the prospects for insights into forward physics at the Large Hadron 
Collider (LHC). We also explain the main electromagnetic processes that gov- 
ern the longitudinal shower evolution. Armed with these two principal shower 
ingredients and motivation from the underlying physics, we describe the dif- 
ferent methods proposed to distinguish primary species. In the last part, we 
outline how ultrahigh energy cosmic ray interactions can be used to probe new 
physics beyond the electroweak scaled 



^Lecture notes from the 6th CERN-Latin- American School of High-Energy Physics, Natal, Brazil, March - April, 2011. 
http : //physicschool . web . cern . ch/PhysicSchool/CLASHEP/CLASHEP201 1/. 
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1 Setting the stage 

For biological reasons our perception of the Universe is based on the observation of photons, most 

trivially by staring at the night-sky with our bare eyes. Conventional astronomy covers many orders of 
magnitude in photon wavelengths, from 10 4 cm radio-waves to 1CT 14 cm gamma rays of GeV energy. 



This 60 octave span in photon frequency allows for a dramatic expansion of our observational capacity 
beyond the approximately one octave perceivable by the human eye. Photons are not, however, the 
only messengers of astrophysical processes; we can also observe cosmic rays and neutrinos (and maybe 
also gravitons in the not-so-distant future). Particle astronomy may be feasible for neutral particles or 
possibly charged particles with energies high enough to render their trajectories magnetically rigid. This 
new astronomy can probe the extreme high energy behavior of distant sources and perhaps provide a 
window into optically opaque regions of the Universe. In these lectures we will focus attention on the 
highest energy particles ever observed. These ultrahigh energy cosmic rays (UHECRs) carry about seven 
orders of magnitude more energy than the LHC beam. We will first discuss the requirements for cosmic 
ray acceleration and propagation in the intergalactic space. After that we will discuss the properties of 
the particle cascades initiated when UHECRs interact in the atmosphere. Finally we outline strategies 
to search for physics beyond the highly successful but conceptually incomplete Standard Model (SM) of 
weak, electromagnetic, and strong interactions. 

Before proceeding, we pause to present our notation. Unless otherwise stated, we work with 
natural (particle physicist's) Heaviside-Lorentz (HL) units with 

h = c = k = £Q = HQ = l. 

The fine structure constant is a = e 2 / (ineohc) ~ 1/137. All SI units can then be expressed in electron 
Volt (eV), namely 

1 m ~ 5.1 x 10 6 eV" 1 , 1 s ~ 1.5 x 10 15 eV" 1 , 1 kg ~ 5.6 x 10 35 eV , 

1 A ~ 1244 eV , 1 G ~ 1.95 x 10~ 2 eV 2 , 1 K ~ 8.62 x 10" 5 eV . 
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Table 1: Symbols and units of most common quantities. 



symbol quantity unit 

J diffuse flux GeV _i cm -2 s _1 sr~ 

<3> diffuse neutrino flux (upper bound) GeV -1 cm -2 s _1 sr~ 

(j) point-source neutrino flux (upper bound) GeV -1 cm -2 s _1 

L source luminosity ergs -1 

C integrated luminosity fb -1 = 10 -39 cm 2 

Q source emissivity GeV -1 s _1 

Q source emissivity per volume GeV -1 cm -3 s -1 



Electromagnetic processes in astrophysical environments are often described in terms of Gauss (G) 
units. For a comparison of formulas used in the literature we note the conversion, (e 2 )nL = 47r(e 2 )c, 
(£> 2 )hl = (£> 2 )g/47t and (£? 2 )hl = {E 2 )q/Ati. To avoid confusion we will present most of the formu- 
las in terms of 'invariant' quantities, i.e. eB, eE and the fine- structure constant a. Distances are generally 
measured in Mpc ~ 3.08 x 10 22 m. Extreme energies are sometimes expressed in EeV = 10 18 eV. The 
symbols and units of most common quantities are summarized in Table 1. 

2 Frontiers of multi-messenger astronomy 
2.1 A century of cosmic ray observations 

In 1912 Victor Hess carried out a series of pioneering balloon flights during which he measured the 
levels of ionizing radiation as high as 5 km above the Earth's surface [1]. His discovery of increased 
radiation at high altitude revealed that we are bombarded by ionizing particles from above. These CR 
particles are now known to consist primarily of protons, helium, carbon, nitrogen and other heavy ions 
up to iron. Measurements of energy and isotropy showed conclusively that one obvious source, the Sun, 
is not the main source. Only below 100 MeV kinetic energy or so, where the solar wind shields protons 
coming from outside the solar system, does the Sun dominate the observed proton flux. Spacecraft 
missions far out into the solar system, well away from the confusing effects of the Earth's atmosphere 
and magnetosphere, confirm that the abundances around 1 GeV are strikingly similar to those found in 
the ordinary material of the solar system. Exceptions are the overabundance of elements like lithium, 
beryllium, and boron, originating from the spallation of heavier nuclei in the interstellar medium. 

Above ~ 10 5 GeV, the rate of CR primaries is less than one particle per square meter per year and 
direct observation in the upper layers of the atmosphere (balloon or aircraft), or even above (spacecraft) 
is inefficient. Only ground-based experiments with large apertures and long exposure times can hope to 
acquire a significant number of events. Such experiments exploit the atmosphere as a giant calorimeter. 
The incident cosmic radiation interacts with the atomic nuclei of air molecules and produces air showers 
which spread out over large areas. Already in 1938, Pierre Auger concluded from the size of the air 
showers that the spectrum extends up to and perhaps beyond 10 6 GeV [2, 3]. In recent years, substantial 
progress has been made in measuring the extraordinarily low flux at the high end of the spectrum. 

There are two primary detection methods that have been successfully used in high exposure exper- 
iments. In the following paragraph we provide a terse overview of these approaches. For an authoritative 
review on experimental techniques and historical perspective see [4]. For a more recent comprehensive 
update including future prospects see [5]. 

The size of an extensive air shower (EAS) at sea-level depends on the primary energy and arrival 
direction. For showers of UHECRs, the cascade is typically several hundreds of meters in diameter 
and contains millions of secondary particles. Secondary electrons and muons produced in the decay of 
pions may be observed in scintillation counters or alternatively by the Cherenkov light emitted in water 
tanks. The separation of these detectors may range from 10 m to 1 km depending on the CR energy 
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Fig. 1: Particles interacting near the top of the atmosphere initiate an electromagnetic and hadronic cascade. Its 
profile is shown on the right. The different detection methods are illustrated. Mirrors collect the Cerenkov and 
nitrogen fluorescent light, arrays of detectors sample the shower reaching the ground, and underground detectors 
identify the muon component of the shower. From Ref . [6] . 

and the optimal cost-efficiency of the detection array. The shower core and hence arrival direction can 
be estimated by the relative arrival time and density of particles in the grid of detectors. The lateral 
particle density of the showers can be used to calibrate the primary energy. Another well-established 
method of detection involves measurement of the shower longitudinal development (number of particles 
versus atmospheric depth, shown schematically in Fig. 1) by sensing the fluorescence light produced via 
interactions of the charged particles in the atmosphere. The emitted light is typically in the 300 - 400 nm 
ultraviolet range to which the atmosphere is quite transparent. Under favorable atmospheric conditions, 
EASs can be detected at distances as large as 20 km, about 2 attenuation lengths in a standard desert 
atmosphere at ground level. However, observations can only be done on clear Moonless nights, resulting 
in an average 10% duty cycle. 

In these lectures we concentrate on results from the High Resolution Fly's Eye experiment and the 
Pierre Auger Observatory to which we will often refer as HiRes and Auger. HiRes [7] was an up-scaled 
version of the pioneer Fly's Eye experiment [8]. The facility was comprised of two air fluorescence 
detector sites separated by 12.6 km. It was located at the US Army Dugway Proving Ground in the 
state of Utah at 40.00° N, 113° W, and atmospheric depth of 870 g/cm 2 . Even though the two detectors 
(HiRes-I and HiRes-II) could trigger and reconstruct events independently, the observatory was designed 
to measure the fluorescence light stereoscopically. The stereo mode allows accurate reconstruction of the 
shower geometry (with precision of 0.4°) and provides valuable information and cross checks about the 
atmospheric conditions at the time the showers were detected. HiRes-I and HiRes-II collected data until 
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April 2006 for an accumulated exposure in stereoscopic mode of 3, 460 hours . The monocular mode 
had better statistical power and covered a much wider energy range. 

The Pierre Auger Observatory [9] is designed to measure the properties of EASs produced by CRs 
at the highest energies, above about 10 9 GeV. It features a large aperture to gather a significant sample 
of these rare events, as well as complementary detection techniques to mitigate some of the systematic 
uncertainties associated with deducing properties of CRs from air shower observables. 

The observatory is located on the vast plain of Pampa Amarilla, near the town of Malargue in 
Mendoza Province, Argentina (35.1° - 35.5° S, 69.6° W, and atmospheric depth of 875 g/cm 2 ). The 
experiment began collecting data in 2004, with construction of the baseline design completed by 2008. 
As of October 2010, Auger had collected in excess of 20, 000 km 2 sr yr in exposure, significantly more 
exposure than other cosmic ray observatories combined. Two types of instruments are employed. Particle 
detectors on the ground sample air shower fronts as they arrive at the Earth's surface, while fluorescence 
telescopes measure the light produced as air shower particles excite atmospheric nitrogen. 

The surface array [10] comprises 1, 600 surface detector (SD) stations, each consisting of a tank 
filled with 12 tons of water and instrumented with 3 photomultiplier tubes which detect the Cherenkov 
light produced as particles traverse the water. The signals from the photomultipliers are read out with 
flash analog to digital converters at 40 MHz and timestamped by a GPS unit, allowing for detailed study 
of the arrival time profile of shower particles. The tanks are arranged on a triangular grid with a 1.5 km 
spacing, covering about 3, 000 km 2 . The surface array operates with close to a 100% duty cycle, and the 
acceptance for events above 10 9 5 GeV is nearly 100%. 

The fluorescence detector (FD) system [11] consists of 4 buildings, each housing 6 telescopes 
which overlook the surface array. Each telescope employs an 11 m 2 segmented mirror to focus the 
fluorescence light entering through a 2.2 m diaphragm onto a camera which pixelizes the image using 
440 photomultiplier tubes. The photomultiplier signals are digitized at 10 MHz, providing a time profile 
of the shower as it develops in the atmosphere. The FD can be operated only when the sky is dark and 
clear, and has a duty cycle of 10-15%. In contrast to the SD acceptance, the acceptance of FD events 
depends strongly on energy [12], extending down to about 10 9 GeV. 

The two detector systems provide complementary information, as the SD measures the lateral 
distribution and time structure of shower particles arriving at the ground, and the FD measures the lon- 
gitudinal development of the shower in the atmosphere. A subset of showers is observed simultaneously 
by the SD and FD. These "hybrid" events are very precisely measured and provide an invaluable cali- 
bration tool. In particular, the FD allows for a roughly colorimetric measurement of the shower energy 
since the amount of fluorescence light generated is proportional to the energy deposited along the shower 
path; in contrast, extracting the shower energy via analysis of particle densities at the ground relies on 
predictions from hadronic interaction models describing physics at energies beyond those accessible to 
current experiments. Hybrid events can therefore be exploited to set a model-independent energy scale 
for the SD array, which in turn has access to a greater data sample than the FD due to the greater live 
time. 

The Extreme Universe Space Observatory (EUSO) is currently being considered by the Japan 
Aerospace Exploration Agency for a possible payload on the Japanese Experimental Module (JEM) of 
the International Space Station [13]. The mission is currently scheduled for 2017 [14]. The launch 
will be provided by the H-II Transfer Vehicle Kounotori. Looking straight down, the UV telescope of 
JEM-EUSO will have 1 km 2 resolution on the Earth's surface, which provides an angular resolution of 
2.5°. The surface area expected to be covered on Earth is about 160, 000 km 2 , with a duty cicle of order 
of 10%. The detector will also operate in a tilted mode that will increase the viewing area by a factor 
of up to 5, but decreasing its resolution. The telescope will be equiped with devices that measure the 
transparency of the atmosphere and the existence of clouds. 

In the remainder of the section, we describe recent results from HiRes and Auger, including the 



5 



10~ 



CO 



CI 

I 



10 4 - 



> 

CD 

o 

r- 



10-* 




10 



13 



Grigorov a 
JACEE v 
MGU v 
TienShan o 
Tibet07 • 
Akeno o 
CASA/MIA eb 
Hegra o 
Flys Eye • 
Agasa ♦ 
HiResl # 
HiRes2 □ 
Auger SD o 
Auger hybrid a 
Kascade ffi 



10 



14 



10 



15 



10 16 10 17 
E [eV] 



10 J 



10 1 



10 



20 



Fig. 2: All-particle spectrum of cosmic rays. From Ref. [15]. 



measurement of the cosmic ray energy spectrum, composition, and searches for anisotropy in the cosmic 
ray arrival directions. 

2.1.1 Energy spectrum 

The CR spectrum spans over roughly 1 1 decades of energy. Continuously running monitoring using so- 
phisticated equipment on high altitude balloons and ingenious installations on the Earth's surface encom- 
pass a plummeting flux that goes down from 10 4 m -2 s _1 at ^ 1 GeV to 10 -2 km -2 yr _1 at ~ 10 11 GeV. 
Its shape is remarkably featureless, with little deviation from a constant power law (J oc i?~ 7 , with 
7^3) across this large energy range. The small changes in the power index, 7, are conveniently visu- 
alized taking the product of the flux with some power of the energy. In this case the spectrum reveals a 
leg-like structure as it is sketched in Fig. 2. The anatomy of this cosmic leg - its changes in slope, mass 
composition or arrival direction - reflects the various aspects of CR propagation, production and source 
distribution. 

A steepening of the spectrum (7 ~ 2.7 — » 3.1) at an energy of about 10 6 5 GeV is known as the 
cosmic ray knee. Composition measurements in cosmic ray observatories indicate that this feature of the 
spectrum is composed of the subsequent fall-off of Galactic nuclear components with maximal energy 
E/Z [16-18]. This scaling with atomic number Z is expected for particle acceleration in a magnetically 
confining environment, which is only effective as long as the particle's Larmor radius is smaller than the 
size of the accelerator. If this interpretation holds, the Galactic contribution in CRs can not extend much 
further than 10 8 GeV, assuming iron {Z = 26) as the heaviest component. However, the observational 
data at these energies is inconclusive and the end-point of Galactic CRs has not been pinned down so far. 
For a survey of spectral features at lower energies see [19]. 



6 



The onset of an extragalactic contribution could be signaled by the so-called second knee - a 
further steepening (7 ~ 3.1 —> 3.2) of the spectrum at about 10 8 7 GeV. Note that extragalactic CRs are 
subject to collisions with the interstellar medium during their propagation over cosmological distances. 
Depending on the initial chemical composition, these particle- specific interactions will be imprinted in 
the spectrum observed on Earth. It has been argued [20, 21] that an extragalactic proton population 
with a simple power-law injection spectrum may reproduce the spectrum above the second knee. In 
these models the flattening (7 ~ 3.2 —> 2.7) in the spectrum at around 10 9 5 GeV - the so-called ankle 
- can be identified as a "dip" from e + e~ pair production together with a "pile-up" of protons from 
pion photoproduction. However, this feature relies on a proton dominance in extragalactic cosmic rays 
since heavier nuclei like oxygen or iron have different energy loss properties in the cosmic microwave 
background (CMB) and mixed compositions, in general will not reproduce the spectral features [22]. A 
cross-over at higher energies is also feasible: above the ankle the Larmor radius of a proton in the galactic 
magnetic field exceeds the size of the Galaxy and it is generally assumed that an extragalactic component 
dominates the spectrum at these energies [23]. Moreover, the galactic-extragalactic transition ought to be 
accompanied by the appearence of spectral features, e.g. two power-law contributions would naturally 
produce a flattening in the spectrum if the harder component dominates at lower energies. Hence, the 
ankle seems to be the natural candidate for this transition. 

Shortly after the CMB was discovered [24], Greisen [25], Zatsepin, and Kuzmin [26] (GZK) 
pointed out that the relic photons make the universe opaque to CRs of sufficiently high energy. This 
occurs, for example, for protons with energies beyond the photopion production threshold, 



where m p (m n ) denotes the proton (pion) mass and c^cmb ~ 10 -3 eV is a typical CMB photon en- 
ergy. After pion production, the proton (or perhaps, instead, a neutron) emerges with at least 50% of 
the incoming energy. This implies that the nucleon energy changes by an e-folding after a propagation 
distance < (a P1 n 7 y^) -1 ~ 15 Mpc. Here, n 7 « 410 cm -3 is the number density of the CMB photons, 
a vl > 0.1 mb is the photopion production cross section, and y n is the average energy fraction (in the lab- 
oratory system) lost by a nucleon per interaction. For heavy nuclei, the giant dipole resonance (GDR) can 
be excited at similar total energies and hence, for example, iron nuclei do not survive fragmentation over 
comparable distances. Additionally, the survival probability for extremely high energy 10 11 GeV) 7- 
rays (propagating on magnetic fields » 10 -11 G) to a distance d, p(> d) « exp[— d/6.6 Mpc], becomes 
less than 10 -4 after traversing a distance of 50 Mpc. All in all, as our horizon shrinks dramatically for 
energies > 10 11 GeV, one would expect a sudden suppression in the energy spectrum if the CR sources 
follow a cosmological distribution. 

At the beginning of summer 2002, in a pioneering paper Bahcall and Waxman [27] noted that 
the energy spectra of CRs reported by the AGASA, the Fly's Eye, the Haverah Park, the HiRes, and 
the Yakutsk collaborations (see Fig. 3) are consistent with the expected GZK suppression at ~ 3. So- 
level according to the Fly's Eye normalization, increasing up to ~ 8a if the selected normalization 
is that of Yakutsk. Five years later, the HiRes Collaboration reported a suppression of the CR flux 
above E = [5.6 ± 0.5(stat) ± 0.9(syst)] x 10 10 GeV, with 5.3a significance [28]. The spectral index 
of the flux steepens from 2.81 ± 0.03 to 5.1 ± 0.7. The discovery of the GZK suppression has been 
confirmed by the Pierre Auger Collaboration, measuring 7 = 2.69 ± 0.2(stat) ± 0.06(syst) and 7 = 
4.2 ± 0.4(stat) ± 0.06(syst) below and above E = 4.0 x 10 10 GeV, respectively (the systematic 
uncertainty in the energy determination is estimated as 22%) [29]. 

Last year, an updated Auger measurement of the energy spectrum was published [30], corre- 
sponding to a surface array exposure of 12, 790 km 2 sr yr. This measurement, combining both hybrid 
and SD-only events, is shown in Fig. 4. The ankle feature and flux suppression are clearly visible. 
A broken power law fit to the spectrum shows that the break corresponding to the ankle is located at 
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Fig. 3: Comparison of flux measurements scaled by E 3 . Only statistical errors are shown. Shown are the data 
of AGASA, Akeno, Auger, Fly's Eye, Haverah Park, HiRes-MIA, HiRes Fly's Eye, MSU, SUGAR, and Yakutsk. 
Yakutsk T500 (trigger 500) refers to the smaller sub-array of the experiment with 500m detector spacing and 
T1000 (trigger 1000) to the array with 1000m detector distance. In case several analyses of the same data set are 
available, only the most recent results are included in the plot. The shaded area, depicting the results of the analysis 
of the Haverah Park data, accounts for some systematic uncertainties by assuming extreme elemental compositions, 
either fully iron or proton dominated. The highest energy point (Fly's Eye monocular observation) corresponds to 
the highest energy event. For sake of clarity, upper limits are not shown. The data of the MSU array are included 
to show the connection of the high-energy measurements to lower energy data covering the knee of the cosmic-ray 
spectrum. FromRef. [19]. 



log 10 (£7/eV) = 18.61 ± 0.01 with 7 = 3.26 ± 0.04 before the break and 7 = 2.59 ± 0.02 after it. The 
break corresponding to the suppression is located at \og 10 (E/eV) = 19.46 ±0.03. Compared to a power 
law extrapolation, the significance of the suppression is greater than 20a. 

2.1.2 Primary species 

Unfortunately, because of the highly indirect method of measurement, extracting precise information 
from EASs has proved to be exceedingly difficult. The most fundamental problem is that the first gen- 
erations of particles in the cascade are subject to large inherent fluctuations and consequently this limits 
the event-by-event energy resolution of the experiments. In addition, the center-of-mass (cm.) energy 
of the first few cascade steps is well beyond any reached in collider experiments. Therefore, one needs 
to rely on hadronic interaction models that attempt to extrapolate, using different mixtures of theory and 
phenomenology, our understanding of particle physics. 

The longitudinal development has a well defined maximum, usually referred to as X max , which 
increases with primary energy as more cascade generations are required to degrade the secondary particle 
energies. Evaluating X max is a fundamental part of many of the composition studies done by detecting 
air showers. For showers of a given total energy E, heavier nuclei have smaller X max because the shower 
is already subdivided into A nucleons when it enters the atmosphere. The average depth of maximum 
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Fig. 4: Combined spectrum from Auger (hybrid and SD events) and the stereo spectrum HiReS. The Auger 
systematic uncertainty of the flux scaled by E 3 , due to the uncertainty of the energy scale of 22%, is indicated by 
arrows. The results of the two experiments are consistent within systematic uncertainties. From Ref. [30]. 

(I max ) scales approximately as \n(E/A) [31]. Therefore, since (X max ) can be determined directly from 
the longitudinal shower profiles measured with a fluorescence detector, the composition can be extracted 
after estimating E from the total fluorescence yield. Indeed, the parameter often measured is Diq, the 
rate of change of (X max ) per decade of energy. 

Photons penetrate quite deeply into the atmosphere due to decreased secondary multiplicities and 
suppression of cross sections by the Landau-Pomeranchuk-Migdal (LPM) effect [32,33]. Indeed, it is 
rather easier to distinguish photons from protons and iron than protons and iron are to distinguish from 
one another. For example, at 10 10 GeV, the (X max ) for a photon is about 1000 g/cm 2 , while for protons 
and iron the numbers are 800 g/cm 2 and 700 g/cm 2 , respectively. 

Searches for photon primaries have been conducted using both the surface and fluorescence istru- 
ments of Auger. While analysis of the fluorescence data exploits the direct view of shower development, 
analysis of data from the surface detector relies on measurement of quantities which are indirectly related 
to the X max , such as the signal risetime at 1000 m from the shower core and the curvature of the shower 
front. Presently, the 95% CL upper limits on the fraction of CR photons above 2, 3, 5, 10, 20, and 
40 x 10 9 GeV are 3.8%, 2.4%, 3.5%, 2.0%, 5.1%, and 31%, respectively. Further details on the analysis 
procedures can be found in [34-36]. In Fig. 5 these upper limits are compared with predictions of the 
cosmogenic photon flux. 

In Fig. 6 we show the variation of (X max ) with primary energy as measured by several experi- 
ments. Interpreting the results of these measurements relies on comparison to the predictions of high 
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Fig. 5: Upper limits on the photon fraction reported by the Pierre Auger Collaboration. A prediction for the 
cosmogenic photon flux is also shown for comparison (details on the calculation are given Sec. 2.5.2) The color 
and line coding is the same as the one in Fig. B.l. This figure is courtesy of Markus Ahlers. 



energy hadronic interaction models. As one can see in Fig. 6, there is considerable variation in pre- 
dictions among the different interaction models. For 1.6 x 10 9 GeV < E < 6.3 x 10 10 GeV, the 
HiRes data are consistent with a constant elongation rate d(X max ) /d(log(E)) = 47.9 ± 6.0(stat) ± 
3.2(sys) g/cm 2 /decade [37]. The inference from the HiRes data is therefore a change in cosmic ray 
composition, from heavy nuclei to protons, at E ~ 10 8 6 GeV [38]. This is an order of magnitude lower 
in energy than the previous crossover deduced from the Fly's Eye data [39]. On the other hand, Auger 
measurements, interpreted with current hadronic interaction models, seem to favor a mixed (protons + 
nuclei) composition at energies above 10 8 6 GeV [40]. However, uncertainties in the extrapolation of the 
proton-air interaction - cross section [41] and elasticity and multiplicity of secondaries [42] - from ac- 
celerator measurements to the high energies characteristic for air showers are large enough to undermine 
any definite conclusion on the chemical composition. 

2.1.3 Distribution of arrival directions 

There exists "lore" that convinces us that the highest energy CRs observed should exhibit trajectories 
which are relatively unperturbed by galactic and intergalactic magnetic fields. Hence, it is natural to 
wonder whether anisotropy begins to emerge at these high energies. Furthermore, if the observed flux 
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Fig. 6: The compilation of fluorescence-based measurements of the mean X max of various experiments is com- 
pared with predictions from EAS simulations using three event generators of hadronic interactions (SIBYLL, 
QGSJET, and DPMJET). The QGSJET predictions on the shower- to- shower fluctuations of the depth of maximum 
are indicated by the shaded (cross-hatched) area for proton (iron) primaries. From Ref. [19]. 



suppression is the GZK effect, there is necessarily some distance, 0(100 Mpc), beyond which cosmic 
rays with energies near 10 11 GeV will not be seen. Since the matter density within about 100 Mpc is 
not isotropic, this compounds the potential for anisotropy to emerge in the UHECR sample. On the one 
hand, if the distribution of arrival directions exhibits a large-scale anisotropy, this could indicate whether 
or not certain classes of sources are associated with large-scale structures (such as the Galactic plane 
or the Galactic halo). On the other hand, if cosmic rays cluster within a small angular region or show 
directional alignment with powerful compact objects, one might be able to associate them with isolated 
sources in the sky. 

CR air shower detectors which experience stable operation over a period of a year or more can have 
a uniform exposure in right ascension, a. A traditional technique to search for large-scale anisotropics is 
then to fit the right ascension distribution of events to a sine wave with period 2ii/m (m th harmonic) to 
determine the components (x, y) of the Rayleigh vector [43] 

2 N 2 N 

x=— r y^ j w i cos(ra OLi) , V = U] i sin(m a { ) , (2) 

i=i i=i 

where the sum runs over the number of N events in the considered energy range, J\f = J2iLi w % * s 
the normalization factor, and the weights, W{ = are the reciprocal of the relative exposure, 

cj, given as a function of the declination, Si [44]. The m th harmonic amplitude of N measurements of 
ai is given by the Rayleigh vector length 1Z = [x 2 + y 2 ) 1 / 2 , and the phase is (/? = arctan(y/x). 
The expected length of such a vector for values randomly sampled from a uniform phase distribution is 
T^o = 2/ \[M . The chance probability of obtaining an amplitude with length larger than that measured is 
p(> H) = e~ k °, where ko = 1Z 2 /1Zq. To give a specific example, a vector of length ko > 6.6 would 
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be required to claim an observation whose probability of arising from random fluctuation was 0.0013 
(a "3cr" result) [45]. For a given CL, upper limits on the amplitude can be derived using a distribution 
drawn from a population characterized by an anisotropy of unknown amplitude and phase 



CL= J- 



[2 1 ds (Tls\ ( s 2 + K 2 /2 \ 

where To is the modified Bessel function of the first kind with order zero and a — \/2/J\[ [43]. 

The first harmonic amplitude of the CR right ascension distribution can be directly related to the 
amplitude a^ of a dipolar distribution of the form J(a, S) = Jo(l + add-u), where u and d respectively 
denote the unit vector in the direction of an arrival direction and in the direction of the dipole. Setting 
m = 1, we can rewrite x, y and J\f as: 

x = — \ dS da cos S J(a, 5) u(S) cos a, 



2 



u min 

pSmax P^TT 

/ dS da cos S J (a, S) u(S) sin a, (4) 

Af = d5 da cos S J (a, S) oj(S) . 

In (4) we have neglected the small dependence on right ascension in the exposure. Next, we write the 
angular dependence in J (a, 5) as d • u — cos 5 cos (5o cos(o: — ao) + sin 5 sin So, where a^ and 5$ are 
the right ascension and declination of the apparent origin of the dipole, respectively. Performing the a 
integration in (4) it follows that 

Ad± 



11 = 

where 



1 + Bd\\ 



(5) 



J d5cj(5) cos 2 5 f d5 (jj(5) cos 5 sin 5 

J d8 uj(S) cos S ' J dS uj(8) cos 6 

d|| = a^sin^o is the component of the dipole along the Earth rotation axis, and d± = a^ cos So is 
the component in the equatorial plane [46]. The coefficients A and B can be estimated from the data 
as the mean values of the cosine and the sine of the event declinations. For example, for the Auger 
data sample we have A = (cos S) ~ 0.78 and B = (sin S) ~ —0.45. For a dipole amplitude a^ the 
measured amplitude of the first harmonic in right ascension 1Z thus depends on the region of the sky 
observed, which is essentially a function of the latitude of the observatory 4ite> and the range of zenith 
angles considered. In the case of a small Bd\\ factor, the dipole component in the equatorial plane d± is 
obtained as d± ~1Z/ (cos S). The phase tp corresponds to the right ascension of the dipole direction a$. 
For a fixed number of arrival directions, the RMS error in the amplitude, Aa^ « 1.57V -1 / 2 , has little 
dependence on the amplitude [44]. 1 

In Fig. 7 we show upper limits and measurements of d± from various experiments together with 
some predictions from UHECR models of both galactic and extragalactic origin. The AGAS A Collabo- 
ration reported a correlation of the CR arrival directions to the Galactic plane at the 4<r level [49]. The 
energy bin width which gives the maximum fco-value corresponds to the region 10 8 9 GeV - 10 9 3 GeV 

AGASA r 

where ko = 11.1, yielding a chance probability of p(> 1Z E _ EeV ) ~ 1.5 x 10 . The recent results 



1 A point worth noting at this juncture: A pure dipole distribution is not possible because the cosmic ray intensity cannot 
be negative in half of the sky. A "pure dipole deviation from isotropy" means a superposition of monopole and dipole, with 
the intensity everywhere > 0. An approximate dipole deviation from isotropy could be caused by a single strong source if 
magnetic diffusion or dispersion distribute the arrival directions over much of the sky. However, a single source would produce 
higher-order moments as well. An example is given in Sec. 2.4. 
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reported by the Pierre Auger Collaboration are inconsistent with those reproted by the AGASA Col- 
laboration [48]. If the galactic/extragalactic transition occurs at the ankle, UHECRs at 10 9 GeV are 
predominantly of galactic origin and their escape from the Galaxy by diffusion and drift motions are 
expected to induce a modulation in this energy range. These predictions depend on the assumed galactic 
magnetic field model as well as on the source distribution and the composition of the UHECRs. Two 
alternative models are displayed in Fig. 7, corresponding to different geometries of the halo magnetic 
fields [50]. The bounds reported by the Pierre Auger Collaboration already exclude the particular model 
with an antisymmetric halo magnetic field (A) and are starting to become sensitive to the predictions of 
the model with a symmetric field (S). The predictions shown in Fig. 7 are based on the assumption of 
predominantly heavy composition in the galactic component [51]. Scenarios in which galactic protons 
dominate at 10 9 GeV would typically predict a larger anisotropy. Alternatively, if the structure of the 
magnetic fields in the halo is such that the turbulent component dominates over the regular one, purely 
diffusion motions may confine light elements of galactic origin up to ~ 10 9 GeV, and may induce an 
ankle-like feature at higher energy due to the longer confinement of heavier elements [52]. Typical sig- 
natures of such a scenario in terms of large scale anisotropics are also shown in Fig. 7 (dotted line). The 
corresponding amplitudes are challenged by the current sensitivity of Auger. On the other hand, if the 
transition is taking place at lower energies, say around the second knee, UHECRs above 10 9 GeV are 
dominantly of extragalactic origin and their large scale distribution could be influenced by the relative 
motion of the observer with respect to the frame of the sources. If the frame in which the UHECR 
distribution is isotropic coincides with the CMB rest frame, a small anisotropy is expected due to the 
Compton-Getting effect [53]. Neglecting the effects of the galactic magnetic field, this anisotropy would 
be a dipolar pattern pointing in the direction ao ~ 168° with an amplitude of about 0.6% [54], close to 
the upper limits set by the Pierre Auger Collaboration. The statistics required to detect an amplitude of 
0.6% at 99% CL is ~ 3 times the published Auger sample [48]. 

The right harmonic analyses are completely blind to intensity variations which depend only on 
declination. Combining anisotropy searches in a over a range of declinations could dilute the results, 
since significant but out of phase Rayleigh vectors from different declination bands can cancel each 
other out. An unambiguous interpretation of anisotropy data requires two ingredients: exposure to the 
full celestial sphere and analysis in terms of both celestial coordinates [44]. 

One way to increase the chance of success in finding out the sources of UHECRs is to check 
for correlations between CR arrival directions and known candidate astrophysical objects. To calculate 
a meaningful statistical significance in such an analysis, it is important to define the search procedure 
a priori in order to ensure it is not inadvertently devised especially to suit the particular data set after 
having studied it. With the aim of avoiding accidental bias on the number of trials performed in selecting 
the cuts, the Auger anisotropy analysis scheme followed a pre-defined process. First an exploratory data 
sample was employed for comparison with various source catalogs and for tests of various cut choices. 
The results of this exploratory period were then used to design prescriptions to be applied to subsequently 
gathered data. 

Based on the results of scanning an exposure of 4, 390 km 2 sr yr, a prescription was designed to test 
the correlation of events having energies E > 5.5 x 10 10 GeV with objects in the Veron-Cetty & Veron 
(VCV) catalog of Active Galactic Nuclei (AGNs). The prescription called for a search of 3.1° windows 
around catalog objects with redshifts z < 0.0018. The significance threshold set in the prescription was 
met in 2007, when the exposure more than doubled and the total number of events reached 27, with 9 
of the 13 events in the post-prescription sample correlating [55,56]. For a sample of 13 events from 
an isotropic distribution, the probability that 9 or more correlate by chance with an object in the AGN 
catalogue (subject to cuts on the exposure weighted fraction of the sky within the opening angles and the 
redshift) is less than 1%. This corresponds to roughly a 2.5<r effect. In the summer of 2008, the HiRes 
Collaboration applied the Auger prescription to their data set and found no significant correlation [57]. 
One has to exercise caution when comparing results of different experiments with potentially different 
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Fig. 7: Upper limits on the anisotropy amplitude of first harmonic as a function of energy from Auger, EAS-TOP, 
AGASA, KASCADE and KASCADE-Grande experiments. Also shown are the predictions up to 1 EeV from two 
different galactic magnetic field models with different symmetries (A and S), the predictions for a purely galactic 
origin of UHECRs up to a few times 10 10 GeV (Gal), and the expectations from the Compton-Getting effect 
for an extragalactic component isotropic in the CMB rest frame (C-G XGal). Below 1 EeV, due to variations 
of the event counting rate arising from atmospheric effects, the Pierre Auger Collaboration adopted the East/West 
method [47], which is two times less efficient but doesn't require correction for trigger efficiency. From Ref. [48] 



energy scales, since the analysis involves placing an energy cut on a steeply falling spectrum. More 
recently, the Pierre Auger Collaboration published an update on the correlation results from an exposure 
of 20, 370 km 2 sr yr (collected over 6 yr but equivalent to 2.9 yr of the nominal exposure/yr of the full 
Auger), which contains 69 events with E > 5.5 x 10 10 GeV [58]. A skymap showing the locations of all 
these events is displayed in Fig. 8. For a physical signal one expects the significance to increase as more 
data are gathered. In this study, however, the significance has not increased. A 3a effect is not neccesarily 
cause for excitement; of every 100 experiments, you expect about one 3a effect. Traditionally, in particle 
physics there is an unwritten 5a rule for "discovery." One should keep in mind though that in the case of 
CR physics we do not have the luxury of controlling the luminosity. 

A number of other interesting observations are described in [58], including comparisons with other 
catalogs as well as a specific search around the region of the nearest active galaxy, Centaurus A (Cen A). 
It is important to keep in mind that these are all a posteriori studies, so one cannot use them to determine 
a confidence level for anisotropy as the number of trials is unknown. A compelling concentration of 
events in the region around the direction of Cen A has been observed. As one can see in Fig. 9, the 
maximum departure from isotropy occurs for a ring of 18° around the object, in which 13 events are 
observed compared to an expectation of 3.2 from isotropy. There are no events coming from less than 
18° around M87, which is almost 5 times more distant than Cen A and lies at the core of the Virgo 
cluster. As shown in Fig. 8 the Auger exposure is 3 times smaller for M87 than for Cen A. Using these 
two rough numbers and assuming equal luminosity, one expects 75 times fewer events from M87 than 
from Cen A. Hence, the lack of events in this region is not completely unexpected. 
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Fig. 8: Correlation of the arrival directions of UHECR with AGNs from the VCV catalog. The shaded part of the 
sky is not visible by Auger. The gray squares are the AGNs within z less than 0.018. The Auger events are shown 
with circles. The first 27 events are half filled. The 13 HiRes events are shown with black dots. The thin lines 
show the six regions of the sky to which Auger has equal exposure. The wide gray line is the supergalactic plane. 
FromRef. [5]. 

The Centaurus cluster lies 45 Mpc behind Cen A. An interesting question then is whether some 
of the events in the 18° circle could come from the Centaurus cluster rather than from Cen A. This does 
not appear likely because the Centaurus cluster is farther away than the Virgo cluster and for compara- 
blel CR luminosities one would expect a small fraction of events coming from Virgo. Furthermore, the 
events emitted by Cen A and deflected by magnetic fields could still register as a correlation due to the 
overdense AGN population lying behind Cen A, resulting in a spurious signal [59]. 




Fig. 9: Cumulative number of events with E > 55 EeV as a function of angular distance from the direction of Cen 
A. The bands correspond to the 68%, 95%, and 99.7% dispersion expected for an isotropic flux. From Ref. [58]. 



15 



In summary, the inaugural years of data taking at the Pierre Auger Observatory have yielded 
a large, high-quality data sample. The enormous area covered by the surface array together with an 
excellent fluorescence system and hybrid detection techniques have provided us with large statistics, 
good energy resolution, and solid control of systematic uncertainties. Presently, Auger is collecting some 
7, 000 km 2 sr yr of exposure each year, and is expected to run for 2 more decades. New detector systems 
are being deployed, which will lower the energy detection threshold down to 10 8 GeV. An experimental 
radio detection program is also co-located with the observatory and shows promising results. As always, 
the development of new analysis techniques is ongoing, and interesting new results can be expected. 

2.2 Origin of ultrahigh energy cosmic rays 

It is most likely that the bulk of the cosmic radiation is a result of some very general magneto-hydrodynamic 
(MHD) phenomenon in space which transfers kinetic or magnetic energy into cosmic ray energy. The 
details of the acceleration process and the maximum attainable energy depend on the particular physical 
situation under consideration. There are basically two types of mechanism that one might invoke. The 
first type assumes the particles are accelerated directly to very high energy (VHE) by an extended electric 
field [60]. This idea can be traced back to the early 1930's when Swann [61] pointed out that betatron 
acceleration may take place in the increasing magnetic field of a sunspot. These so-called "one-shot" 
mechanisms have been worked out in greatest detail, and the electric field in question is now generally 
associated with the rapid rotation of small, highly magnetized objects such as neutron stars (pulsars) or 
AGNs. Electric field acceleration has the advantage of being fast, but suffers from the circumstance that 
the acceleration occurs in astrophysical sites of very high energy density, where new opportunities for 
energy loss exist. Moreover, it is usually not obvious how to obtain the observed power law spectrum in a 
natural way, and so this kind of mechanism is not widely favored these days. The second type of process 
is often referred to as statistical acceleration, because particles gain energy gradually by numerous en- 
counters with moving magnetized plasmas. These kinds of models were mostly pioneered by Fermi [62]. 
In this case the E~ 2 spectrum emerges very convincingly. However, the process of acceleration is slow, 
and it is hard to keep the particles confined within the Fermi engine. In this section we first provide a 
summary of statistical acceleration based on the simplified version given in Ref. [63]. For a more detailed 
and rigorous discussion, the reader is referred to [64]. After reviewing statistical acceleration, we turn to 
the issue of the maximum achievable energy within diffuse shock acceleration and explore the viability 
of some proposed UHECR sources. 

2.2.1 Fermi acceleration at shock waves 

In his original analysis, Fermi [62] considered the scattering of CRs on moving magnetized clouds. 
The right panel of Fig. 10 shows a sketch of these encounters. Consider a CR entering into a single 
cloud with energy E{ and incident angle 0i with the cloud's direction undergoing diffuse scattering on 
the irregularities in the magnetic field. After diffusing inside the cloud, the particle's average motion 
coincides with that of the gas cloud. The energy gain by the particle, which emerges at an angle Of 
with energy Ef, can be obtained by applying Lorentz transformations between the laboratory frame 
(unprimed) and the cloud frame (primed). In the rest frame of the moving cloud, the CR particle has a 
total initial energy 

E'i = Idoud /^cioud cos 0i) , (6) 

where r c i ou d and /3 c ioud = Kioud/c are the Lorentz factor and velocity of the cloud in units of the speed 
of light, respectively. In the frame of the cloud we expect no change in energy (E[ = E'j), because 
all the scatterings inside the cloud are due only to motion in the magnetic field (so-called collisionless 
scattering). There is elastic scattering between the CR and the cloud as a whole, which is much more 
massive than the CR. Transforming to the laboratory frame we find that the energy of the particle after 
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its encounter with the cloud is 



E f = r cloud E) (1 + /3 cloud cos 6 f ) . (7) 
The fractional energy change in the laboratory frame is then 

AE _ E f - Ei 1 - /3 d oud cos 0i + /3 doud cos Of - /^ loud cos 0; cos Of 

Inside the cloud the CR direction becomes randomized and so (cos Of) = 0. The collisionless scattered 
particle will gain energy in a head-on collision (0^ > tt/2) and lose energy by tail-end (0{ < tt/2) 
scattering. The net increase of its energy is a statistical effect. The average value of cos 0{ depends on the 
relative velocity between the cloud and the particle. The probability P per unit solid angle of having a 
collision at angle 0{ is proportional to (v — Kioud cos 0$), where v is the CR speed. In the ultrarelativistic 
limit, i.e., v ~ c (as seen in the laboratory frame), 

dP 

— oc (1 - /3 doud cos 0^) , (9) 

so 

<cos^) = -^p>. (10) 



Now, inserting Eq. (10) into Eq.(8), one obtains for /3 c i ou d <^ 1> 

3c2loud/3 - 1 - 4 ^ loud . (id 



(AE) _ 1 + /3 c 2 loud /3 4 

— — -L r *-~ / — 

E 1 - ^loud 3 



Note that (AE)/E oc /3 doud , so even though the average magnetic field may vanish, there can still 
be a net transfer of the macroscopic kinetic energy from the moving cloud to the particle. However, 
the average energy gain is very small, because /3 do ud ^ acceleration process is very similar 

to a thermodynamical system of two gases, which tries to come into thermal equilibrium [66]. Corre- 
spondingly, the spectrum of CRs should follow a thermal spectrum which might be in conflict with the 
observed power-law. 

A more efficient acceleration may occur in the vicinity of plasma shocks occurring in astrophysical 
environments [67,68]. Suppose that a strong (nonrelativistic) shock wave propagates through the plasma 
as sketched in the left panel of Fig. 10. Then, in the rest frame of the shock the conservation relations 
imply that the upstream velocity ^ up (ahead of the shock) is much higher than the downstream velocity 
^down (behind the shock). The compression ratio r = ^ U p/^down = ^down/^up can be determined by 
requiring continuity of particle number, momentum, and energy across the shock; here n up (nd 0W n) is 
the particle density of the upstream (downstream) plasma. For an ideal gas the compression ratio can be 
related to the specific heat ratio and the Mach number of the shock. For highly supersonic shocks, r = 
4 [64]. Therefore, in the primed frame stationary with respect to the shock, the upstream flow approaches 
with speed u up = f3 up c = 4/3c/3 and the downstream flow recedes with speed ^ doW n = AiownC = /3c/ 3. 
When measured in the stationary upstream frame, the quantity u = u np — ^ doW n = /3c is the speed of 
the shocked fluid and u up = /3 s hock is the speed of the shock. Hence, because of the converging flow 
- whichever side of the shock you are on, if you are moving with the plasma, the plasma on the other 
side of the shock is approaching you with velocity u - to first order there are only head-on collisions for 
particles crossing the shock front. The acceleration process, although stochastic, always leads to a gain in 
energy. In order to work out the energy gain per shock crossing, we can visualize magnetic irregularities 
on either side of the shock as clouds of magnetized plasma of Fermi's original theory. By considering 
the rate at which CRs cross the shock from downstream to upstream, and upstream to downstream, one 
finds (cos Oi) = —2/3 and (cos Of) = 2/3. Hence, Eq. (8) can be generalized to 

( AE ) ~ip = 4 ^np-^down j (12) 

E 3 3 c 
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Note this is first order in = u/c, and is therefore more efficient than Fermi's original mechanism. 

An attractive feature of Fermi acceleration is its prediction of a power-law flux of CRs. Consider 
a test-particle with momentum p in the rest frame of the upstream fluid (see Fig. 10). The particle's 
momentum distribution is isotropic in the fluid rest frame. For pitch angles tt/2 < Oi < tt relative to 
the shock velocity vector (see Fig. 10) the particle enters the downstream region and has on average 
the relative momentum p[l + 2(/3 up — /3down)/3]. Subsequent diffusion in the downstream region 're- 
isotropizes' the particle's momentum distribution in the fluid rest frame. As the particle diffuses back 
into the upstream region (for pitch angles < Of < tt/2) it has gained an average momentum of 
(Ap) jp ~ 4(/3 up — /3down)/3. This means that the momentum gain of a particle per time is proportional 
to its momentum, 

P = P Again • (13) 

On the other hand, the loss of particles from the acceleration region is proportional to their number, 

N = -N/t loss . (14) 
Therefore, taking the ratio (13)/(14) we first obtain 

dN/dp= -aN/p, (15) 

and after integration N(p) oc p~ a , with a = t ga { n /t\ oss . If the acceleration cycle across the shock 
takes the time At we have already identified At/t ga i n = {Ap) /p ~ 4(/3 up — /?down)/3. For the loss 
of relativistic particles one finds At/t\ oss ~ 4/3 up . Therefore, a ~ 3/3 up /(/3 up — /3d 0W n) = 3r/(r — 1), 
yielding a ~ 4 for highly supersonic shocks and a > 4 otherwise. The energy spectrum N(E) oc 
is related to the momentum spectrum by dEN(E) = Aiip 2 &pN(p) and hence 7 = a — 2 > 2. The 
steeply falling spectrum of CRs with 7^3 seems to disfavor supersonic plasma shocks. However, for 
the comparison of these injection spectra with the flux of CRs observed on Earth, one has to consider 
particle interactions in the source and in the interstellar medium. This can have a great impact on the 
shape, as we will discuss in Sees. 2.4 and 2.5.2. 

In general, the maximum attainable energy of Fermi's mechanism is determined by the time scale 
over which particles are able to interact with the plasma. For the efficiency of a "cosmic cyclotron" 
particles have to be confined in the accelerator by its magnetic field B over a sufficiently long time scale 
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compared to the characteristic cycle time. The Larmor radius of a particle with charge Ze increases with 
its energy E according to 



1 E 



Ana ZeB 

-¥Gi*)(£)"V 

The particle's energy is limited as its Larmor radius approaches the characteristic radial size i? S ource of 
the source 

E -^ z (^)(^) xio9Gev - <iv ' 

This limitation in energy is conveniently visualized by the 'Hillas plot' [60] shown in Fig. 1 1 where the 
characteristic magnetic field B of candidate cosmic accelerators is plotted against their characteristic size 
R. It is important to stress that in some cases the acceleration region itself only exists for a limited period 
of time; for example, supernovae shock waves dissipate after about 10 4 yr. In such a case, Eq. (17) would 
have to be modifed accordingly. Otherwise, if the plasma disturbances persist for much longer periods, 
the maximum energy may be limited by an increased likelihood of escape from the region. A look at 
Fig. 11 reveals that the number of sources for the extremely high energy CRs around 10 12 GeV is very 
sparse. For protons, only radio galaxy lobes and clusters of galaxies seem to be plausible candidates. 
For nuclei, terminal shocks of galactic superwinds originating in the metally-rich starburst galaxies are 
potential sources [69]. Exceptions may occur for sources which move relativistically in the host-galaxy 
frame, in particular jets from AGNs and gamma-ray bursts (GRBs). In this case the maximal energy 
might be increased due to a Doppler boost by a factor ^ 30 or ^ 1000, respectively. 

For an extensive discussion on the potential CR-emitting- sources shown in Fig. 11, see e.g. [70]. 
Two of the most attractive examples are discussed next. 



2.2.2 AGNs 

AGNs are composed of an accretion disk around a central super-massive black hole and are sometimes 
associated with jets terminating in lobes which can be detected in radio. One can classify these objects 
into two categories: radio-quiet AGN with no prominent radio emission or jets and radio-loud objects 
presenting jets. 

Fanaroff-Riley II (FRII) galaxies [71] are the largest known dissipative objects (non-thermal sources) 
in the cosmos. Localized regions of intense synchrotron emission, known as "hot spots," are observed 
within their lobes. These regions are presumably produced when the bulk kinetic energy of the jets 
ejected by a central active nucleus is reconverted into relativistic particles and turbulent fields at a "work- 
ing surface" in the head of the jets [72]. Specifically, the speed ^head ~ ^jet [1 + (^e/^jet) 1 ^ 2 ] -1 , with 
which the head of a jet advances into the intergalactic medium of particle density n e can be obtained 
by balancing the momentum flux in the jet against the momentum flux of the surrounding medium; 
where nj e t and Uj e t are the particle density and the velocity of the jet flow, respectively (for relativistic 
corrections, see [73]). For n e > nj e t, ^jet > ^head so that that the jet decelerates. The result is the 
formation of a strong collisionless shock, which is responsible for particle reacceleration and magnetic 
field amplification [74]. The acceleration of particles up to ultrarelativistic energies in the hot spots is the 
result of repeated scattering back and forth across the shock front, similar to that discussed in Sec. 2.2.1. 
The particle deflection in this mechanism is dominated by the turbulent magnetic field with wavelength 
k equal to the Larmor radius of the particle concerned [75]. A self-consistent (although possibly not 
unique) specification of the turbulence is to assume that the energy density per unit of wave number of 
MHD turbulence is of the Kolmogorov type, I(k) oc A: -5 / 3 , just as for hydrodynamical turbulence [76]. 
With this in mind, to order of magnitude accuracy using effective quantities averaged over upstream (jet) 
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Fig. 11: The "Hillas plot" for various CR source candidates (blue areas). Also shown are jet-frame parameters 
for blazers, gamma-ray bursts, and microquasars (purple areas). The corresponding point for the LHC beam is 
also shown. The red dashed lines show the lower limit for accelerators of protons at the CR knee (~ 10 6 5 GeV), 
CR ankle (~ 10 9,5 GeV) and the GZK suppression (~ 10 10,6 GeV). The dotted gray line is the upper limit 
from synchrotron losses and proton interactions in the cosmic photon background (R ^> 1 Mpc). The grey area 
corresponds to astrophysical environments with extremely large magnetic field energy that would be gravitationally 
unstable. FromRef. [65]. 



and downstream (hot spot) conditions (considering that downstream counts a fraction of 4/5 [75]) the 
acceleration timescale at a shock front is found to be [77] 
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is the Kolmogorov diffusion coefficient, U is the ratio of turbulent to ambient magnetic energy density 
in the region of the shock (of radius R), and B is the total magnetic field strength. 



20 



The subtleties surrounding the conversion of particles kinetic energy into radiation provide ample 
material for discussion. The most popular mechanism to date relates 7-ray emission to the development 
of electromagnetic cascades triggered by secondary photomeson products that cool instanstaneously via 
synchrotron radiation. The characteristic single photon energy in synchrotron radiation emitted by an 
electron is 

^ n = (£) 1/2 £^f§ ~ 5 - 4 x 10-2 ^° ^/ Eey ) 2 Tev • ™ 

For a proton this number is (m p /m e ) 3 ~6x 10 9 times smaller. 

The acceleration process will then be efficient as long as the energy losses by synchrotron radiation 
and/or photon-proton interactions do not become dominant. The synchrotron loss time for protons is 
given by [78] 

Ts y n ~ 2TD2 ' ( 21 ) 

where or and T = E/(m p c 2 ) are the Thomson cross section and Lorentz factor, respectively. Con- 
sidering an average cross section a 1P [79] for the three dominant pion-producing interactions, jp —> 
p7T°, 7p —> n7r + , 7p —> p7T + 7T~ , the time scale of the energy losses, including synchrotron and photon 
interaction losses, reads [80] 

6tt mi C 3 - T SV n 

E' 1 = syn A , (22) 



loss ~ a T m 2 B 2 (1 + Aa) 1 + Aa 

where a stands for the ratio of photon to magnetic energy densities and A gives a measure of the relative 
strength of ^p interactions versus the synchrotron emission. Note that the second channel involves the 
creation of ultrarelativistic neutrons that can readily escape the system. For typical hot spot conditions, 
the number density of photons per unit energy interval follows a power-law spectrum 

AGN/ a J (N /CJ ) (CJ/CJO)- 2 U Q <UJ<UJ* 

^ W-\ otherwise ^ 

where No is the normalization constant and uo and uo* correspond to radio and gamma rays energies, 
respectively. The ratio of photon to magnetic energy density is then 

a ~ B 2 /8tt (M) 
and A is only weakly dependent on the properties of the source 

A= a JL (m p /m e ) 2 ^^n> lQx 1q5 ^ 20Q (25) 
ax m(o;*/a;o) 

The maximum attainable energy can be obtained by balancing the energy gains and losses [81] 

E 20 = 1.4 x 10 5 B-W f% (1 + Aa)^ , (26) 

where E = 10 20 E20 eV and R = R\^ pc 1 kpc. It is of interest to apply the acceleration conditions to the 
nearest AGN. 

At only 3.4 Mpc distance, Cen A is a complex FRI radio-loud source identified at optical frequen- 
cies with the galaxy NGC 5128 [82]. Radio observations at different wavelengths have revealed a rather 
complex morphology shown in Fig. 12. It comprises a compact core, a jet (with subluminal proper mo- 
tions /3jet ~ 0.5 [85]) also visible at X-ray frequencies, a weak counterjet, two inner lobes, a kpc-scale 
middle lobe, and two giant outer lobes. The jet would be responsible for the formation of the northern 
inner and middle lobes when interacting with the interstellar and intergalactic media, respectively. There 
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Fig. 12: Optical image of Cen A (UK 48-inch Schmidt) overlaid with radio contours (black, VLA [83]). H.E.S.S. 
VHE best fit position with la statistical errors (blue cross) and VHE extension 95% CL upper limit (white dashed 
circle) are also shown. From Ref. [84]. 

appears to be a compact structure in the northern lobe, at the extrapolated end of the jet. This structure 
resembles the hot spots such as those existing at the extremities of FRII galaxies. However, at Cen A it 
lies at the side of the lobe rather than at the most distant northern edge, and the brightness contrast (hot 
spot to lobe) is not as extreme [86]. Estimates of the radio spectral index of synchrotron emission in 
the hot spot and the observed degree of linear polarization in the same region suggests that the ratio of 
turbulent to ambient magnetic energy density in the region of the shock is U ~ 0.4 [87]. The broadband 
radio- to-X-ray jet emission yields an equipartition magnetic field B^q ~ 100 [88]. 2 The radio- visible 
size of the hot spot can be directly measured from the large scale map R\^ pc ~ 2 [89]. The actual size 
can be larger because of uncertainties in the angular projection of this region along the line of sight. 3 Re- 
placing these fiducial values in (17) and (26) we conclude that if the ratio of photon to magnetic energy 
density a < 0.4, it is plausible that Cen A can accelerate protons up E « 2 x 10 11 GeV. 

EGRET observations of the gamma ray flux for energies > 100 MeV allow an estimate L 7 ~ 
10 41 erg s _1 for Cen A [90]. This value of L 7 is consistent with an earlier observation of photons in the 
TeV-range during a period of elevated X-ray activity [91], and is considerably smaller than the estimated 
bolometric luminosity Lbol ~ 10 43 erg s _1 [82]. Recent data from H.E.S.S. have confirmed Cen A as a 
TeV 7-ray emitting source [84]. Extrapolating the spectrum measured with EGRET in the GeV regime 

2 The usual way to estimate the magnetic field strength in a radio source is to minimize its total energy. The condition of 
minimum energy is obtained when the contributions of the magnetic field and the relativistic particles are approximately equal 
(equipartition condition). The corresponding B-field is commonly referred to as the equipartion magnetic field. 

3 For example, an explanation of the apparent absence of a counterjet in Cen A via relativistic beaming suggests that the 
angle of the visible jet axis with respect to the line of sight is at most 36° [86], which could lead to a doubling of the hot spot 
radius. It should be remarked that for a distance of 3.4 Mpc, the extent of the entire source has a reasonable size even with this 
small angle. 
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to VHEs roughly matches the H.E.S.S. spectrum, though the softer end of the error range on the EGRET 
spectral index is preferred. More recent data from Fermi-LAT established that a large fraction (> 1/2) 
of the total > 100 MeV emission from Cen A emanates from the lobes [92]. For values of B in the 
100 fiG range, substantial proton synchrotron cooling is suppressed, allowing production of high energy 
electrons through photomeson processes. The average energy of synchrotron photons scales as [93] 

{E^ n ) -0.29£; yn , (27) 

and therefore, to account for the observed TeV photons Cen A should harbor a population of ultra- 
relativistic electrons, with E e ~ 10 9 GeV. We further note that this would require the presence of protons 
with energies between one and two orders of magnitude larger, since the electrons are produced as 
secondaries. 



2.2.3 GRBs 

GRBs are flashes of high energy radiation that can be brighter, during their brief existence, than any other 
source in the sky. The bursts present an amazing variety of temporal profiles, spectra, and timescales [94]. 
Our insights into this phenomenon have been increased dramatically by BATSE observations of over 
2000 GRBs, and more recently, by data from SWIFT. 

There are several classes of bursts, from single-peaked events, including the fast rise and expo- 
nential decaying (FREDs) and their inverse (anti-FREDs), to chaotic structures [95]. There are well 
separated episodes of emission, as well as bursts with extremely complex profiles. Most of the bursts are 
time asymmetric, but some are symmetric. Burst timescales range from about 30 ms to several minutes. 

The GRB angular distribution appears to be isotropic, suggesting a cosmological origin [96]. Fur- 
thermore, the detection of "afterglows" — delayed low energy (radio to X-ray) emission — from GRBs 
has confirmed this via the redshift determination of several GRB host-galaxies [97]. 

The 7-ray luminosity implied by cosmological distances is astonishing: L 1 ~ 10 52 erg/s. The 
most popular interpretation of the GRB -phenomenology is that the observable effects are due to the 
dissipation of the kinetic energy of a relativistic expanding plasma wind, a "fireball" [98]. Although the 
primal cause of these events is not fully understood, it is generally believed to be associated with the 
core collapse of massive stars (in the case of long duration GRBs) and stellar collapse induced through 
accretion or a merger (short duration GRBs) [99]. 

The very short timescale observed in the light curves indicates an extreme compactness (i.e. dis- 
tance scale comparable to a light-ms: ro ~ 10 7 cm) that implies a source which is initially opaque 
(because of 77 pair creation) to 7-rays 

r 77 ~ r n« RB a T ~ ~ 10 15 , (28) 

1 1 j 47T ro c e 7 

where n^ RB is number density of photons at the source and e 7 ~ 1 MeV is the characteristic photon 
energy. 

The high optical depth creates the fireball: a thermal plasma of photons, electrons, and positrons. 
The radiation pressure on the optically thick source drives relativistic expansion (over a time scale ro / c), 
converting internal energy into the kinetic energy of the inflating shell. As the source expands, the optical 
depth is reduced. If the source expands with a Lorentz factor T, the energy of photons in the source frame 
is smaller by a factor T compared to that in the observer frame, and most photons may therefore be below 
the pair production threshold. Baryonic pollution in this expanding flow can trap the radiation until most 
of the initial energy has gone into bulk motion with Lorentz factors of T ~ 10 2 — 10 3 [100]. The 
kinetic energy can be partially converted into heat when the shell collides with the interstellar medium 
or when shocks within the expanding source collide with one another. The randomized energy can then 
be radiated by synchrotron radiation and inverse Compton scattering yielding non-thermal bursts with 
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timescales of seconds at large radii, r > 10 12 cm, beyond the Thompson sphere. Charged particles may 
be efficiently accelerated to ultrahigh energies in the fireball's internal shocks, hence GRBs are often 
considered as potential sources of UHECRs [101]. 

Coburn and Boggs [102] reported the detection of polarization, a particular orientation of the 
electric-field vector, in the 7-rays observed from a burst. The radiation released through synchrotron 
emission is highly polarized, unlike in other previously suggested mechanisms such as thermal emission 
or energy loss by relativistic electrons in intense radiation fields. Thus, polarization in the 7-rays from a 
burst provides direct evidence in support of synchrotron emission as the mechanism of 7-ray production 
(see also [103]). 

Following Hillas criterion, the Larmor radius tl should be smaller than the largest scale Zgrb 
over which the magnetic field fluctuates, since otherwise Fermi acceleration will not be efficient. One 
may estimate Zgrb as follows. The comoving time, i.e. the time measured in the fireball rest frame, is 
t = r/Tc. Hence, the plasma wind properties fluctuate over comoving scale length up to /grb ~ r/T, 
because regions separated by a comoving distance larger than r/T are causally disconnected. Moreover, 
the internal energy is decreasing because of the expansion and thus it is available for proton acceleration 
(as well as for 7-ray production) only over a comoving time t. The typical acceleration time scale is 
then [101] 

^ B ~ T ± • (29) 

Equation (29) sets a constraint on the required comoving magnetic field strength, and the Larmor radius 
tl — E' /eB — E/FeB, where E' — E/Y is the proton energy measured in the fireball frame. This 
constraint sets a lower limit to the magnetic field carried by the wind, which may be expressed as 

<B > Q M T h^k^ (3Q) 
Ce £52 

where T = 10 2 5 r2.5, L 1 = 10 52 Ls2 erg s _1 . Here, (b is the fraction of the wind energy density which 
is carried by the magnetic field, Anr 2 T 2 (B 2 / '8tt) = (bL, and ( e is the fraction of wind energy carried 
by shock accelerated electrons. Note that because the electron energy is lost radiatively, L 1 « ( e L. 

The dominant energy loss process in this case is synchrotron cooling. Therefore, the condition 
that the synchrotron loss time of Eq. (21) be smaller than the acceleration time sets the upper limit on the 
magnetic field strength [101] 

" ° w m T 2.5 ^20 

Since the equipartition field is inversely proportional to the radius r, this condition may be satisfied 
simultaneously with (30) provided that the dissipation radius is large enough, i.e. 

r > 10 12 T^ 2 E% cm. (32) 

The high energy protons also lose energy in interaction with the wind photons (mainly through pion 
production). It can be shown, however, that this energy loss is less important than the synchrotron 
energy loss [101]. 

In summary, a dissipative ultra-relativistic wind, with luminosity and variability time implied by 
GRB observations, satisfies the constraints necessary to accelerate protons to energy > 10 11 GeV, pro- 
vided that T > 100, and the magnetic field is close to equipartition with electrons. 



B < 3 x 10 5 r| 5 E~ 2 G. (31) 



2.3 Energy losses of baryonic cosmic rays on the pervasive radiation fields 
2.3.1 Opacity of the CMB to UHECR protons 

Ultrahigh energy protons degrade their energy through Bethe-Heitler (BH) pair production (pj —> 
pe + e~) and photopion production (pj — >► ttN), each successively dominating as the proton energy in- 
creases. The fractional energy loss due to interactions with the cosmic background radiation at a redshift 
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z = is determined by the integral of the nucleon energy loss per collision multiplied by the probability 
per unit time for a nucleon collision in an isotropic gas of photons [104]. This integral can be explicitly 
written as follows, 

IdE c ^ [ u ™ 1 . . nJu) 

~Elt = 2^^ J * ^ i r/2T du ~W ' (33) 

where uo r is the photon energy in the rest frame of the nucleon, and yj is the inelasticity, i.e. the average 
fraction of the energy lost by the photon to the nucleon in the laboratory frame for the jth reaction 
channel. (Here the laboratory frame is the one in which the CMB is at a temperature ^ 2.7 K.) The sum 
is carried out over all channels, n 7 (uj)du stands for the number density of photons with energy between 
uj and duj, aj(uj r ) is the total cross section of the jth interaction channel, T is the usual Lorentz factor of 
the nucleon, and uj m is the maximum energy of the photon in the photon gas. 

Pair production and photopion production processes are only of importance for interactions with 
the 2.7 K blackbody background radiation [105]. Collisions with optical and infrared photons give a 
negligible contribution. Therefore, for interactions with a blackbody field of temperature T, the photon 
density is that of a Planck spectrum, so the fractional energy loss is given by 

where uo j is the threshold energy for the jth reaction in the rest frame of the nucleon. 

At energies E <C m e m p /kT = 2.1 x 10 9 GeV (i.e., uj r /m e - 2 C 1), when the reaction takes 
place on the photons from the high energy tail of the Planck distribution, the fraction of energy lost in 
one collision and the cross section can be approximated by the threshold values 

V m = 2 ^ , (35) 

I Up 

and 



W e ~ 2 ) ' (36) 

where a is the fine structure constant and ro is the classical radius of the electron [105]. The fractional 
energy loss due to pair production for E < 10 9 GeV is then, 

1 (dE\ lQcm e 2 ( kT\ 3 (FkT\ 2 / m e \ 

~e \ii) Bn = — m P ar ° \m) UtJ ex n-fkf)- (37) 

At higher energies (E > 10 10 GeV) the characteristic time for the energy loss due to pair production 
is t ~ 5 x 10 9 yr [106]. In this energy regime, the photopion reactions —> pn and pj —> n+n on 
the tail of the Planck distribution give the main contribution to proton energy loss. The cross sections of 
these reactions are well known and the kinematics is simple. 

Photopion production turns on at a photon energy in the proton rest frame of 145 MeV with a 
strongly increasing cross section at the A (1232) resonance, which decays into the one pion channels 
7T + n and 7T°p. With increasing energy, heavier baryon resonances occur and the proton might reappear 
only after successive decays of resonances. The most important channel of this kind is p^ A ++ 7r~ 
with intermediate A ++ states leading finally to A ++ — >► p?r + . A ++ examples in this category are the 
A(1620) and A(1700) resonances. The cross section in this region can be described by either a sum or 
a product of Breit-Wigner distributions over the main resonances produced in 7 collisions considering 
tvN, tvtvN and KA (A Ntv) final states [107]. At high energies, 3.0 GeV < u r < 183 GeV, the 
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E\ = vi^A_ ^ (40) 



CERN-HERA and COMPAS Groups have made a fit to the p^ cross section [108]. The parameterization 
is 

a n (u r ) = A + B In 2 (^) + C In (^) mb , (38) 

where A = 0.147 ± 0.001, B = 0.0022 ± 0.0001, and C = -0.0170 ± 0.0007. In this energy range, 
the cr to tai(^7) is to a good approximation identical to cr to tai(P7)- 

We turn now to the kinematics of photon-nucleon interactions. The inelasticity y n depends not 
only on the outgoing particles but also on the kinematics of the final state. Nevertheless, averaging over 
final state kinematics leads to a good approximation of y n . The cm. system quantities (denoted by 
*) are determined from the relativistic invariance of the square of the total 4-momentum p^ of the 
photon-proton system. This invariance leads to the relation 

s = (cj* + E*) 2 = m 2 p + 2m p u r . (39) 

The cm. system energies of the particles are uniquely determined by conservation of energy and mo- 
mentum. For reactions mediated by resonances one can assume a decay, which in the cm. frame is 
symmetric in the forward and backward directions with respect to the collision axis (given by the in- 
coming particles). For instance, we consider single pion production via the reaction p^ —> A —> pn. 
Here, 

(s + m\ - m 2 ) 
2^~s 

Thus, the mean energy of the outgoing proton is 

/ F *finah = (s + m 2 A - ml) (m| + m 2 - ml) 

[ p ' 2^m A 2m A ' K } 

/pfinalv E (g ~ ml + mj) {m\ - m 2 + m 2 ) 
{E ? ) = 7 2^ 2^ • (42) 

The mean inelasticity y n = 1 — ( (i? final ) /E) of a reaction that provides a proton after n resonance decays 
can be obtained by straightforward generalization of Eq. (42), and is given by 

(43) 

where m# denotes the mass of the i tn resonant system of the decay chain, mM the mass of the associ- 
ated meson, trr = is the total energy of the reaction in the cm., and rriR n the mass of the nucleon. 
For multi-pion production the case is much more complicated because of the non-trivial final state kine- 
matics. However, it is well established experimentally [109] that, at very high energies (y/s > 3 GeV), 
the incoming particles lose only one-half their energy via pion photoproduction independently of the 
number of pions produced, y n ~ 1/2. This is the "leading particle effect". 

For < 2 GeV, the best maximum likelihood fit to Eq. (34) with the exponential behavior 



or in the lab frame 
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E l¥ 



Aexp[-B/E], (44) 



derived from the values of cross section and fractional energy loss at threshold, gives [110] 

A = (3.66 ± 0.08) x 10" 8 yr -1 , B = (2.87 ± 0.03) x 10 11 GeV . (45) 

The fractional energy loss due to production of multipion final states at higher cm. energies (y/s > 
3 GeV) is roughly a constant, 

) =C= (2.42 ± 0.03) x 10~ 8 yr" 1 . (46) 
E \dt J n 
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Fig. 13: Energy attenuation length of protons in the intergalactic medium. Note that after a distance of ~ 100 Mpc, 
or propagation time ~ 3 x 10 8 yr, the mean energy is essentially independent of the initial energy of the protons, 
with a critical energy around 10 11 GeV. From Ref. [110]. 



From the values determined for the fractional energy loss, it is straightforward to compute the energy 
degradation of UHECRs in terms of their flight time. This is given by, 

At - Ei (B/E) + Ei (B/Eq) = , for 10 10 GeV < E < 10 12 GeV , (47) 

and 

E(t) = E exp[- C t ] , for E > 10 12 GeV , (48) 

where Ei is the exponential integral [111]. Figure 13 shows the proton energy degradation as a function 
of the mean flight distance. Notice that, independent of the initial energy of the nucleon, the mean energy 
values approach 10 11 GeV after a distance of « 100 Mpc. 

2.3.2 Photonuclear interactions 

The relevant mechanisms for the energy loss that extremely high energy nuclei suffer during their trip 
to Earth are: Compton interactions, pair production in the field of the nucleus, photodisintegration, and 
hadron photoproduction. The Compton interactions have no threshold energy. In the nucleus rest-frame, 
pair production has a threshold at ^ 1 MeV, photodisintegration is particularly important at the peak of 
the GDR (15 to 25 MeV), and photomeson production has a threshold energy of ~ 145 MeV. 

Compton interactions result in only a negligibly small energy loss for the nucleus given by [112] 
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where p 1 is the energy density of the ambient photon field in eV cm -3 , E is the total energy of the 
nucleus in eV, and Z and A are the atomic number and weight of the nucleus. The energy loss rate 
due to photopair production is Z 2 /A times higher than for a proton of the same Lorentz factor [113], 
whereas the energy loss rate due to photomeson production remains roughly the same. The latter is 
true because the cross section for photomeson production by nuclei is proportional to the mass number 
A [114], while the inelasticity is proportional to 1/A However, it is photodisintegration rather than 
photopair and photomeson production that determines the energetics of UHECR nuclei. During this 
process some fragments of the nuclei are released, mostly single neutrons and protons. Experimental 
data of photonuclear interactions are consistent with a two-step process: photoabsorption by the nucleus 
to form a compound state, followed by a statistical decay process involving the emission of one or more 
nucleons. 

Following the conventions of Eq. (33), the disintegration rate with production of i nucleons is 
given by [115] 

1 f°° n (uj) f 2Vuj 

R Ai=^2 J du 1 J1 J d0J r U r (7 A i{Ur) (50) 

with cjAi the cross section for the interaction. 

The photoabsorption cross section roughly obeys the Thomas-Reiche-Kuhn (TRK) dipole sum 

rule 

f°° N Z 

E d = / a{uj r )duj r = 59.8 —— MeVmb , (51) 
Jo A 

where N = A — Z is the number of neutrons. (Indeed, this integral is experimentally ^ 20 — 30% larger, 

e.g. for 56 Fe, 1, 020 mb-MeV for the left hand side, 22% larger than the right hand side [116].) These 

cross sections contain essentially two regimes. At uo r < 30 MeV there is the domain of the GDR where 

disintegration proceeds mainly by the emission of one or two nucleons. A Gaussian distribution in this 

energy range is found to adequately fit the cross section data [112]. At higher energies, the cross section 

is dominated by multinucleon emission and is approximately flat up to u r ~ 150 MeV. Specifically, 

- 2) O(30 - ay) e -2K-eo0 2 /A? /iSd e(a; r - 30) 

° M = WA~i + 120 ' (52) 

for i=\, 2, and 

/iS d e(a; r -30) .... 

°M = ' (53) 

for i > 2 [112]. Here, W is a normalization factor given by 

W = g) 2 [$(^2(30 - e 0l )/A z ) + $(^2(6o, - 2)/A<)] , 

<3>(x) is the error function, and Q(x) the Heaviside step function. The dependence of the width A^, the 
peak energy eoi, the branching ratio fi, and the dimensionless integrated cross section ^ are given in 
Ref. [112] for isotopes up to 56 Fe. 

The photon background relevant for nucleus disintegration consists essentially of photons of the 
2.7 K CMB. The background of optical radiation turns out to be of (almost) no relevance for UHECR 
propagation. The cosmic infrared background (CIB) radiation [117] 

^p- = 1.1 x ID" 4 f^-)" 2 ' 5 cm-SeV" 1 , (54) 
duo VeV/ 



only leads to sizeable effects far below 10 11 GeV and for time-scales O (10 17 s) [118]. 

By substituting Eqs. (52) and (53) into Eq. (50) the photodisintegration rates on the CMB can be 
expressed as integrals of two basic forms. The first one is 

*i5/r 



2r 2 7T 2 ft 3 C 2 



ri-o/L roc 
Jl/Y J15/T 



(55) 
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where the functions J and J' are given by the expressions, 



J = 



+ 



eoiAi \ <S>(V2{2Tuj - eoi)/A<) + <S>(V2(e 0i - 2)/A< 



(56) 



A. 



-2((e *-2)/A,) 2 _ -2((2ru;-e i)/Ai 



'}• 



and 



J' = J-eoiAi [$(^2(30 - eoi)/A<) + *( V2(e 0i - 2)/A< 

2 



(57) 



+ 



9< 



e -2((e ,-2)/A,) 2 _ e -2((30-e 0l )/A, 



The second basic integral is of the form 



I 2 = (tt^VJ-Va 



ii 5 /r e"/ fcT - 1 ~ V ~r) L/r 



e uj/kT i 



(58) 



With this in mind, Eq. (50) can be re- written as [1 19] 



^J f kT\n(l-e- 15 / rkT ) 



-±e-*&M (£) VZ f <S 3 + | JCkT [HI ~ e~^ kT ) - ln(l - e~V™) 



TT \ 1/2 A? 



A 



+ 



120 



r 2 5 4 + 15 2 fcrin(l-e- 15 / rfer ) 



(59) 



with A, Si, and K as given in Table 2. Summing over all the possible channels for a given number of 
nucleons, one obtains the effective nucleon loss rate Ra = iR-Ai- The effective nucleon loss rate for 
light elements, as well as for those in the carbon, silicon and iron groups can be scaled as [1 12] 



dA 
~dt 



dA 
~dt 



56 Fe 



#56 



A 
56 



with the photodisintegration rate (59) parametrized by [120] 

R 56 (T) = 3.25 x 10- 6 r- a643 exp(-2.15 x 10 10 /r) s" 1 
forT G [1.0 x 10 9 ,36.8 x 10 9 ], and 

i? 5 e(r) = i.59 x io _i2 r- 

forT G [3.68 x 10 10 , 10.0 x 10 10 ]. 



-12-T-0.0698 -1 



(60) 



(61) 



(62) 



EXERCISE 2.1 Approximating the cross section in Eq. (52) by the single pole of the Narrow- 
Width Approximation [121] 

ga{oj') = ir&o r ° ( f R 8(u>' - luq) , (63) 



show that for interactions with the CMB photons 



Ra 



4r 2 vr 



In (l - e-<l 2TT ) 



(64) 
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Table 2: Series and functions of Eq. (59). 



A 




Si 


E7 =1 kTf 1 exp[B z ] + 15V8/A0 - #(£ + >/8/A<)} 


s 2 


E7 =1 ^i" 1 exp{- J yrfcT}[$( v / 2(2 - eoiJ/Ai) - $(^(30 - eoO/AO] 


S3 


1 exp[^] + 15v/8/Ai) - $(B + VS/Ai)} 


S 4 


El, exp{-15j/rfcT}[(fcT/j)(15/r) 2 + (fcT/j) 2 (15/r) + {kT/jf] 


B 


jAi/FkTVS2 - 2e 0i /V2Ai 


K 


e 0i Ai <S>(V2(e 0i - 2)/A;) + (A;/2) 2 exp{-2(£ 0i - 2) 2 /A 2 } 



where a /A = 1.45 x lCT 27 cm 2 , r GDR = 8 MeV, and e' = 42.65A~ 21 (0.925A 2 433 ) MeV, for 
A > 4 (A < 4) [122]. Verify that for 56 Fe this solution agrees to within 20% with the parametrization 
given in Eq. (61). 

For photodisintegration, the averaged fractional energy loss results equal the fractional loss in 
mass number of the nucleus, because the nucleon emission is isotropic in the rest frame of the nucleus. 
During the photodisintegration process the Lorentz factor of the nucleus is conserved, unlike the cases 
of pair production and photomeson production processes which involve the creation of new particles that 
carry off energy. The total fractional energy loss is then 

_ 1 dE _ 1 dY R 

~E^tt~f^tt + A' ( } 

For u r < 145 MeV the reduction in T comes from the nuclear energy loss due to pair production. The 7- 
ray momentum absorbed by the nucleus during the formation of the excited compound nuclear state that 
precedes nucleon emission is (9(10 -2 ) times the energy loss by nucleon emission [123]. For T > 10 10 
the energy loss due to photopair production is negligible, and thus 



E{t) - 938 A(t) T MeV 



Eq exp 



56 



(66) 



Figure 14 shows the energy of the heaviest surviving nuclear fragment as a function of the prop- 
agation time, for initial iron nuclei. The solid curves are obtained using Eq. (66), whereas the dashed 
and dotted-dashed curves are obtained by means of Monte Carlo simulations [118]. One can see that 
nuclei with Lorentz factors above 10 10 cannot survive for more than 10 Mpc. For these distances, the 
approximation given in Eq. (66) always lies in the region which includes 95% of the Monte Carlo simula- 
tions. When the nucleus is emitted with a Lorentz factor Tq < 5 x 10 9 , pair production losses start to be 
relevant, significantly reducing the value of T as the nucleus propagates distances of (9(100 Mpc). The 
effect has a maximum for Tq « 4 x 10 9 but becomes small again for Tq < 10 9 , for which appreciable 
effects only appear for cosmological distances (> 1000 Mpc), see for instance [118]. 

Note that Eq. (66) imposes a strong constraint on the location of nucleus-sources: less than 1% 
of iron nuclei (or any surviving fragment of their spallations) can survive more than 3 x 10 14 s with an 
energy > 10 1L5 GeV For straight line propagation, this represents a maximum distance of ^ 3 Mpc. 



2.4 Diffuse propagation of protons in a magnetized Local Supercluster 

In addition to the interactions with the radiation fields permeating the universe, baryonic CRs suffer 
deflection and delay in magnetic fields, effects which can camouflage their origins. For example, the 
regular component of the Galactic magnetic field can distort the angular images of CR sources: the flux 
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Fig. 14: The energy of the surviving fragment (To = 4 x 10 9 , To = 2 x 10 10 ) vs. propagation time obtained 
using Eq. (66) is indicated with a solid line. Also included is the energy attenuation length obtained from Monte 
Carlo simulations with (dashed) and without (dotted-dashed) pair creation production, for comparison. The region 
between the two dotted lines includes 95% of the simulations. This gives a clear idea of the range of values which 
can result from fluctuations from the average behaviour. It is important to keep in mind that a light propagation 
distance of 1.03 x 10 14 s corresponds to 1 Mpc. From Ref. [124]. 

may appear dispersed around the source or globally translated in the sky with rather small dispersion, viz. 
the deflection for CRs of charge Ze and energy E should not exceed - 10° Z(4x 10 10 GeV/E) [125]. 

One interesting possibility to explain the observed near-isotropy of arrival directions is to envis- 
age a large scale extragalactic magnetic field that can provide sufficient bending to the CR trajectories. 
Surprisingly little is actually known about the extragalactic magnetic field strength. There are some 
measurements of diffuse radio emission from the bridge area between the Coma and Abell superclus- 
ters [126], which under assumptions of equipartition allows an estimate of (9(0.2 — 0.6) /xG for the 
magnetic field in this region. Fields of O(fiG) are also indicated in a more extensive study of 16 low 
redshift clusters [127]. It is assumed that the observed £> -fields result from the amplification of much 
weaker seed fields. However, the nature of the initial week seed fields is largely unknown. There are two 
broad classes of models for seed fields: cosmological models, in which the seed fields are produced in 
the early universe, and astrophysical models, in which the seed fields are generated by motions of the 
plasma in (proto)galaxies. Of particular interest here is the second class of models. If most galaxies lived 
through an active phase in their history, magnetized outflows from their jets and winds would efficiently 
pollute the extragalactic medium. The resulting £>-field is expected to be randomly oriented within cells 
of sizes below the mean separation between galaxies, A# < 1 Mpc. 

Extremely weak unamplified extragalactic magnetic fields have escaped detection up to now. Mea- 
surements of the Faraday rotation in the linearly polarized radio emission from distant quasars [128] 
and/or distortions of the spectrum and polarization properties in the CMB [129, 130] imply upper limits 
on the extragalactic magnetic field strength as a function of the reversal scale. It is important to stress that 
Faraday rotation measurements (RM) sample extragalactic magnetic fields of any origin (out to quasar 
distances), while the CMB analyses set limits only on primordial magnetic fields. The RM bounds de- 
pend significantly on assumptions about the electron density profile as a function of the redshift. When 
electron densities follow that of the Lyman-a forest, the average magnitude of the magnetic field receives 
an upper limit of B ~ 10 -9 G for reversals on the scale of the horizon, and B ~ 10 -8 G for reversal 
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scales on the order of 1 Mpc [131]. As a statistical average over the sky, an all pervading extragalactic 
magnetic field is constrained to be [132] 

B<3x 1(T 7 (Q^/Om)- 1 (Zi/0.72) (A^/Mpc) 1 / 2 G , (67) 

where Q^h 2 ~ 0.02 is the baryon density and h ~ 0.72 is the present day normalized Hubble expansion 
rate. This is a conservative bound because has contributions from neutrons and only electrons in 
ionized gas are relevant to Faraday rotation. 

In the spirit of [133, 134], very recently we proposed that neutron emission from Cen A could 
dominate the observed CR flux above the GZK suppression [135]. Neutrons that are able to decay 
generate proton diffusion fronts in the intergalactic turbulent magnetic plasma. In our calculations we 
assume a strongly turbulent magnetic field: B = 50 nG, Xb ~ 1 Mpc, and largest turbulent eddy 
£ ~ 2tt\b [136]. For energies above the GZK supression, Xb < tl < £ and so the diffusion coefficient 
is given by the Bohm formula [75] 



The evolution of the proton spectrum is driven by the so-called "energy loss-diffusion equation" 

dn(E,r,t) _ d[b(E)n(E , r , t)} 



(68) 



dt dE 



+ V[D(E, r, t) Vn(E, r, *)] + Q(E, t) S 3 (x) , (69) 



where n(E, r, 0) = No S 3 (x). Here, b(E) = dE/dt is the mean rate at which particles lose energy 
and Q(E,t) is the number of protons per unit energy and per unit time generated by the source. For 
the situation at hand, D(E,r,t) = D{E) and hence the second term becomes D(E)V 2 n(E,r,t). 
Idealizing the emission to be uniform with a rate dNo /dt = N tot / r, we have 

Q(E,t) = ^ [Q(t - ton) - 0(t - toff)], (70) 
r 

where / Q(E, r', t') d 3 x' dt' = N tot , O is the Heaviside step function, and t on (t g) is the time since the 
engine turned on (off) its CR production, £ Q ff — ^on = t. For the energy region of interest, the expected 
time delay of the diffuse protons, raday ~ d 2 /D(E), is significantly smaller than the characteristic time 
scale for photopion production derived in Eq. (44). 

If the energy loss term is neglected, the solution to Eq. (69) reads, 

n(E, r,t) = j dt' J d 3 x' G(r -r',t- t') Q(E, r', t'), (71) 

where 

G(r-r',t-t') = [47rD(t-t')]" 3/2 ©(* - O exp{-(r - r') 2 /4D(t - t')} , (72) 

is the Green function [137]. The density of protons at the present time t of energy E at a distance r from 
Cen A, which is assumed to be continuously emitting at a constant spectral rate dNo / dE dt from time 
i on until the present, is found to be [134] 

dn(E, r, t) _ dN 1 [ l ^ e -r 2 /4£>(t-f ) 



dE dE dt [4ttL>(£)] 3 / 2 J tm {t - t'fl 2 

r 2 d 

)r J vi v 3 



dN 1 r dv _ 1/v 

dEdt 4n 3 / 2 D(E)r J V1 v 3 ^ 
dNo 1 



dE dt 4irD(E)r 



I(x) , (73) 
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where we have used the change of variables 



U =4W^)=l (74) 



with x = 4DT on /r 2 , T on = t — t on , and 



du 



1 f°° 

I{x) = -f= / e~ u . (75) 

For T on — >► oc, the density approaches its time-independent equilibrium value n eq . 

As a result of this diffusion the J oc E~ A behavior of the observed CR spectrum reflects a 
dNo/dEdt oc E~ 3 injection in the region of the source cutoff. For S = 3000 km 2 detector like Auger, 
the neutron rate is 



f 



dt And? J El dEdt 

where \{E) ~ 9.2 x 10" 3 E^eV Mpc is the neutron decay length. For the energy interval between E\ = 
55 EeV and E2 = 150 EeV, we calculate the normalization factor using the observation of 2 neutrons 
in 3 yr of the nominal exposure/yr of Auger. We then use this normalization factor to calculate the 
luminosity of the source in the above energy interval. We find L^ E2) = 0.86 x 10 40 erg/s [135]. Next, 
we assume continuity of the spectrum at E\ as it flattens to E~ 2 . Taking the lower bound on the energy 
to be Eq = 10 EeV, we can then fix the luminosity for this interval and find L^ ,El ^ = 2.3 x 10 40 erg/s. 

Adding these, we find the (quasi) bolometric luminosity to be Lq^'^ = 3.2 x 10 40 erg/s, which is about 
a factor of 3 smaller than the observed luminosity in E > 100 MeV 7-rays L 1 « 10 41 erg/s [90]. To 
further constrain the parameters of the model, we evaluate the energy-weighted approximately isotropic 
proton flux at 70 EeV. If the source actively emitted UHECRs for at least 70 Myr, from Eq. (73) we 
obtain 

{E J{E)) = ^)HD{E) dEdt I{x) 

w 1.6 x 10 57 eV 3 km" 2 yr" 1 sr" 1 , (77) 

in agreement with observations [29]. If we assume circular pixel sizes with 3° radii, the neutrons will be 
collected in a pixel representing a solid angle Afi ^ 8.6 x 10" 3 sr. The event rate of (diffuse) protons 
coming from the direction of Cen A is found to be 

- Q AO /^ 2 /Z?4 7 \ ^ - O HQ nr^ + o/.rr. (J 



P 2 dE 
S An / (E 4 J) —- = 0.08 events/yr . 



d t j El 

All in all, in the next 9 yr of operation we expect about 6 direct neutron events against an almost negligible 
background. Note that our model also predicts no excess in the direction of M87, in agreement with 
observations (see Fig. 8). 

We turn now to the discussion of anisotropy. The number of particles with velocity c hitting a unit 
area in a unit time in a uniform gas of density n(E, r , t) is n(E, r , t) c. Due to the gradient in the number 
density with radial distance from the source, the downward flux at Earth per steradian as a function of 
the angle 9 to the source is [133] 

J(0, r, t) = n(i Y' t)c (l + a d cos 6) , (79) 

where 

a d cos6=— — . (80) 

n(E,r,t) c 
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The asymmetry parameter can be found by computing the incoming current flux density j = D^V^n 
as viewed by an observer on Earth, where r = R + r' . We obtain 

n ^0 1 f* dt' -(B?+2R.r'+r' 2 )/4D(t-t') ~( 2R i + 2x 'i) 

dEdt (47rD) 3 / 2 J ton (t - t'f/ 2 4D(i - f) 



(Ri + x't) r dt' 



f 

J ton 



e -rV[4D(t-t')] (gl) 



2(4ttD)3/2 J ton (t-t'f' 2 
Near = 0, using the change of variables in Eq. (74), we obtain 

tj 2 d£dt (47TD) 3 / 2 4D ^ r 2 y A/x 

^ ^° 1 /'(,), (82) 



2vr dEdt r 3 



where 



If 00 

l'(x) = -= / du^e~ u . 
V 71 " Jl/x 



(83) 



l/a 

Finally, taking R x = R x = and R z = r cos # we obtain [134] 



2g(g) /'(*) 

For T on oo, < a < 1. Taking our fiducial values E = 70 EeV, 5 nG = 50, and T on = 70, we find 
ad = 0.29. This is within ~ la of the anisotropy amplitude ad = d±/ cos«5o = 0.25 ± 0.18 obtained 
from the 69 arrival directions, assuming a dipole function for a source model with a maximum value at 
CenA: (a ,6 ) = (201.4°, -43.0°). 

One caveat is that we assumed that neutrons completely dominate the ultrahigh energy Cen A 
emission spectrum; that is 

^ = (JVJ + N$)E- S , with JVJ/JVJ«1. (85) 

This reduces the number of free parameters in the model. The actual proton-to-neuton fraction depends 
on the properties of the source, dominantly by the ratio of photon-to-magnetic energy density. The rela- 
tion (85) results from a ~ 0.4, which implies that photopion production - and not proton leakege from 
the accelerator region - determines the shape of the cutoff in the spectrum at the source. 



EXERCISE 2.2 The assymetry parameter given in Eq. (84) accurately describes the diffuse prop- 
agation of CRs if Nq/Nq <C 1. Show that asymmetry parameter from the neutron decay shell is given 
by 

2D(E) I"(x) 



OLd 



cr I{x) 



where 



(1+k) 2 u 



(86) 



(87) 



and k = \{E)/r. Show that in spite of the complicated nature of Eq.(87), the results for ay are very 
similar to the ones for the primary diffusion front given above. 
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Fig. 15: The nominal arrival directions (a = right ascension, 5 = declination) of SUGAR events with energies 
above 4 x 10 10 GeV [138]. Also shown in solid lines are contour maps indicating the circular areas of the celestial 
sphere centered at Cen A (indicated by *) with 10° and 25° radii. The dashed lines surrounding several of the 
events indicate the uncertainty of each arrival direction; this was found to be about 3° sec 0, where is the shower 
zenith angle. From Ref. [81]. 



An observation demanding of some attention is the distribution of arrival directions of the events 
collected by the SUGAR experiment, shown in Fig. 15. We take the data at face value keeping in 
mind that the techniques available at the time the experiment was conducted were not as refined as the 
techniques currently at our disposal. As one can see in Fig. 15, while the Auger and SUGAR distrbutions 
are not dissimilar, the statistics are rather limited and do not support a firm conclusion. 

As an attentive reader should know by now, the observed excess in the direction of Cen A can 
also be explained using proton directional signals and diffusion of heavy nuclei in a £>-field of about 
1 nG [139]. However, if this were the case, a larger anisotropy (as yet unobserved) should be present 
at EeV energies [140]. Regarding the preceeding discussion, one should note that the 3° window does 
not have an underlying theoretical motivation. Recall that this angular range resulted from a scan of 
parameters to maximize signal significance. Cen A covers an eliptical region spanning about 10° along 
the major axis; see Fig. 12. Therefore, some care is required to select the region of the sky which is most 
likely to maximize the signal-to-noise. 

In summary, existing data is consistent with the hypothesis that Cen A dominates the CR sky at 
the high end of the spectrum [133]. Auger is in a gifted position to explore Cen A and would provide in 
the next 9 yr of operation sufficient statistics to test this hypothesis. The potential detection of neutrons 
at Auger can subsequently be validated by the larger aperture of the JEM-EUSO mission. 

2.5 Ultrahigh energy cosmic neutrinos 
2.5.1 Waxman-Bahcall bound 

It is helpful to envision the CR engines as machines where protons are accelerated and (possibly) per- 
manently confined by the magnetic fields of the acceleration region. The production of neutrons and 
pions and subsequent decay produces neutrinos, 7-rays, and CRs. If the neutrino-emitting source also 
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produces high energy CRs, then pion production must be the principal agent for the high energy cutoff 
on the proton spectrum. Conversely, since the protons must undergo sufficient acceleration, inelastic 
pion production needs to be small below the cutoff energy; consequently, the plasma must be optically 
thin. Since the interaction time for protons is greatly increased over that of neutrons due to magnetic 
confinement, the neutrons escape before interacting, and on decay give rise to the observed CR flux. The 
foregoing can be summarized as three conditions on the characteristic nucleon interaction time scale 7i nt ; 
the neutron decay lifetime r n \ the characteristic cycle time of confinement r cyc i e ; and the total proton 
confinement time r conf : (i) r int > r cyc i e ; (ii) r n > r cyc i e ; (Hi) r int < r conf . The first condition 
ensures that the protons attain sufficient energy. Conditions (i) and (ii) allow the neutrons to escape the 
source before decaying. Condition (Hi) permits sufficient interaction to produce neutrons and neutrinos. 
These three conditions together define an optically thin source [141]. A desirable property to reproduce 
the almost structureless energy spectrum is that a single type of source will produce cosmic rays with a 
smooth spectrum across a wide range of energy. 

The cosmic ray flux above the ankle is often summarized as "one 3 x 10 10 GeV particle per 
kilometer square per year per steradian." This can be translated into an energy flux [142] 

E(EJ CR ) - 3Xl0l ° GeV 



(10 10 cm 2 )(3 x 10 7 s)sr 
= 10" 7 GeVcm- 2 s- 1 sr- 1 . (88) 



From this we can derive the energy density ecR in UHECRs using flux = velocity x density, or 

4tt j dE{EJ CR } = ce C R. (89) 



This leads to 



r-£ max 10 -7 QeV „_ 1Q TeV 



4lT f m ' it; .. rv 

e C R = — / ~^dE ~ 10" 19 ^3 , (90) 

c . h cm z s cm- 3 



taking the extreme energies of the accelerator(s) to be E m [ n — 10 10 GeV and £ ma x = 10 12 GeV. 
The power required for a population of sources to generate this energy density over the Hubble time 

(T H « 10 10 yr) is: e§R°' 10l2] - 5 x 10 44 TeV Mpc~ 3 yr" 1 - 3 x 10 37 erg Mpc~ 3 s" 1 . This works out 
to roughly (i) L « 3 x 10 39 erg s _1 per galaxy, (ii) L « 3 x 10 42 erg s _1 per cluster of galaxies, (Hi) 
L « 2 x 10 44 erg s _1 per active galaxy, or (iv) J Ldt « 10 52 erg per cosmological GRB [142]. The 
coincidence between these numbers and the observed output in electromagnetic energy of these sources 
explains why they have emerged as the leading candidates for the CR accelerators. 

The energy production rate of protons derived professionally, assuming a cosmological distribu- 
tion of sources (with injection spectrum typical of shock acceleration dNo/dE oc E~ 2 ) is [143] 

^CR ' 1012 ' ~ 5 x 1044 er g M PC" 3 yr" 1 . (91) 

This is within a factor of our back-of-the-envelope estimate (1 TeV =1.6 erg). The energy-dependent 
generation rate of CRs is therefore given by 



E 2 d ™ 



.[10 10 ,10 12 ] 



dE ln(10 12 /10 10 ) 

w 10 44 ergMpc" 3 yr- 1 . (92) 

The energy density of neutrinos produced through interactions of these protons can be directly 
tied to the injection rate of CRs 

o dn v 3 ^ o dh 
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where is the fraction of the energy which is injected in protons lost into photopion interactions. The 
factor of 3/8 comes from the fact that, close to threshold, roughly half the pions produced are neutral, 
thus not generating neutrinos, and one quarter of the energy of charged pion decays goes to electrons 
rather than neutrinos. Namely, resonant p~f interactions produce twice as many neutral pions as charged 
pions. Direct pion production via virtual meson exchange contributes only about 20% to the total cross 
section, but is almost exclusively into 7r + . Hence, p~f interactions produce roughly equal numbers of 7r + 
and 7T°. The average neutrino energy from the direct pion decay is {E y ^Y = (1 — r) E n /2 ~ 0.22 E n 
and that of the muon is (E^) 71 = (1 + r) ~ 0.78 E^, where r is the ratio of muon to the pion 

mass squared. In muon decay, since the has about 1/3 of its parent energy, the average muon neutrino 
energy is (E^Y = (1 + r)E n /6 = 0.26 E n . 

The "Waxman-Bahcall bound" is defined by the condition = 1 

« 2.3 x lO^e^&GeVcm^s^sr" 1 , (94) 

where the parameter £ z accounts for the effects of source evolution with redshift, and is expected to 
be ^ 3 [144]. For interactions with the ambient gas (i.e. pp rather than p^ collisions), the average 
fraction of the total pion energy carried by charged pions is about 2/3, compared to 1/2 in the photopion 
channel. In this case, the upper bound given in Eq. (94) is enhanced by 33%. Electron antineutrinos can 
also be produced through neutron /3-decay. This contribution, however, turns out to be negligible (see 
Appendix A for details). 

The actual value of the neutrino flux depends on what fraction of the proton energy is converted to 
charged pions (which then decay to neutrinos), i.e. is the ratio of charged pion energy to the emerging 
nucleon energy at the source. For resonant photoproduction, the inelasticity is kinematically determined 
by requiring equal boosts for the decay products of the A + , giving = E^+ /E n « 0.28, where 
and E n are the emerging charged pion and neutron energies, respectively. For pp —> NN + pions, where 
N indicates a final state nucleon, the inelasticity is « 0.6 [145]. This then implies that the energy carried 
away by charged pions is about equal to the emerging nucleon energy, yielding (with our definition) 
« 1. 

At production, if all muons decay, the neutrino flux consists of equal fractions of v e , and v^. 
Originally, the Waxman-Bahcall bound was presented for the sum of and (neglecting u e ), moti- 
vated by the fact that only muon neutrinos are detectable as track events in neutrino telescopes. Since 
oscillations in the neutrino sector mix the different species, we chose instead to discuss the sum of all 
neutrino flavors. When the effects of oscillations are accounted for, nearly equal numbers of the three 
neutrino flavors are expected at Earth [146]. 



EXERCISE 2.3 The assumption that GRBs are the sources of the observed UHECRs generates 
a calculable flux of neutrinos produced when the protons interact with the fireball photons [147]. In 
the observer's frame, the spectral photon density (GeV -1 cm -3 ) can be adequately parametrized by a 
broken power-law spectrum n^ RB (e 7 ) oc e^, where (3 ~ 1, 2 respectively at energies below and above 

e break _ j Me y [148] Show ^ 

/ inbreak \ — 1 

«b(^ > ^ break ) - ID" 13 [Y^GeV ) Cm " 2 S_1 Sr_1 ' (95) 

where E^ ieak ~5x 10 5 r| 5 (e^ reak /MeV)- 1 GeV. Recall that in our convention e 7 = Te' T where 
is the photon energy measured in the fireball frame. Convince yourself that the non-observation of 
extraterrestrial neutrinos from sources other than the Sun and SN 1987a puts the GRB model of UHECR 
acceleration on probation [149]. 
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If the injected cosmic rays include nuclei heavier than protons, then the neutrino flux expected 
from the cosmic ray sources may be modified. Nuclei undergoing acceleration can produce pions, just 
as protons do, through interactions with the ambient gas, so the Waxman-Bahcall argument would be 
unchanged in this case. However, if interactions with radiation fields dominate over interactions with 
matter, the neutrino flux would be suppressed if the cosmic rays are heavy nuclei. This is because the 
photodisintegration of nuclei dominates over pion production at all but the very highest energies. Defin- 
ing n as the fraction of nuclei heavier than protons in the observed cosmic ray spectrum, the resulting 
neutrino flux is then given by E% $^ a11 « (1 — k) E% [150]. The most up-to-date calculation of k 

combines a double-fit analysis of the energy and elongation rate measurements to constrain the spectrum 
and chemical composition of UHECRs at their sources [151]. Injection models with a wide range of 
chemical composition are found to be consistent with observations. In particular, the data are consistent 
with a proton dominated spectrum with only a small (1 - 10%) admixture of heavy nuclei. 4 

By duplicating the Waxman-Bahcall calculation for Cen A, we obtain an upper limit on the inten- 
sity of neutrinos from the direction of the nearest active galaxy, 

« 5.0 x 10" 9 GeV cm" 2 s" 1 , (96) 

with E v < 10 9 - 5 GeV. For the model introduced in Sec. 2.4, we have 671- ~ 0.28 and therefore a 
prediction for the all-flavor neutrino flux 

El (/)" c f nA « 1.5 x 10~ 9 GeV cm- 2 s" 1 , (97) 

in agreement with the results of [152]. In addition to neutrinos, one also expects a similar flux of photons 
at the source, which also carries unperturbed directional information. However, photons at these energies 
are difficult to dig out from the hughe proton background. 

Although there are no other nearby FRI of this magnitude which can potentially be detected as 
point sources, one can integrate over the estimated FRI population out to the horizon to obtain an upper 
limit for the diffuse FRI neutrino flux. This quantity is given by [153] 



1_ _ . 3 

47T 

1.5 x 10~ 8 GeV cm- 2 s" 1 sr" 1 , (98) 



where TZh — 1 horizon ~ 3 Gpc and n FRI ^8x 10 4 Gpc -3 is the number density [154]. Note that 
this flux is about a factor of 3 smaller than the Waxman-Bahcall upper limit. Hence, the reduction in 
luminosity of the ensemble of neutrino sources roughly compensates for the presence of distant optically 
thin sources whose CR components are hidden by extragalactic magnetic fields. 

The diffuse neutrino flux has an additional component originating in the energy losses of UHECRs 
en route to Earth [155]. The accumulation of these neutrinos over cosmological time is known as the cos- 
mogenic neutrino flux. The GZK reaction chain generating cosmogenic neutrinos is well known [156]. 
The intermediate state of the reaction P7cmb —> n7T + /pn is dominated by the A + resonance, because 
the neutron decay length is smaller than the nucleon mean free path on the CMB. Gamma-rays, pro- 
duced via 7T° decay, subsequently cascade electromagnetically on intergalactic radiation fields through 
e + e~ pair production followed by inverse Compton scattering. The net result is a pile up of 7-rays at 
GeV-TeV energies, just below the threshold for further pair production on the diffuse optical background. 
Meanwhile each tt + decays to 3 neutrinos and a positron; the e + readily loses its energy through inverse 
Compton scattering on the diffuse radio background or through synchrotron radiation in intergalactic 
magnetic fields. As we have seen, the neutrinos carry away about 3/4 of the 7r + energy, therefore the 



4 It is important to stress that the essential results of the analysis in Ref. [151] are not altered by the new Auger data [40]. 
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energy in cosmogenic neutrinos is about 3/4 of that produced in 7-rays. The intensity of 7-ray pile up 
currently provides the most stringent bound on the flux of cosmogenic neutrinos. It is this that we now 
turn to study. 



2.5.2 Boltzmann equation, universal cosmic ray spectrum, and cosmogenic neutrinos 

For a spatially homogeneous distribution of cosmic sources, emitting UHE particles of type i, the comov- 
ing number density, Y{(z, E) = rti(z, E)/(l + z) 3 , is governed by a set of one-dimensional (Boltzman) 
continuity equations, 



Yi = d E (HEYi) + d E {biYi) - V t Y % + ^ J dE 3 ljt Y 3 + Q % , 



(99) 



together with the Friedman-Lemaitre equations describing the cosmic expansion rate H(z) as a function 
of the redshift z. 5 The first term on the r.h.s. describes adiabatic energy losses [6]. The second term 
describes interactions on the cosmic photon backgrounds which can be approximated by continuous 
energy losses. The third and fourth terms describe more general interactions involving particle losses 
(i -» anything) with interaction rate 1^, and particle generation of the form j -» i. The last term on the 
r.h.s., Qi, corresponds to the luminosity density per comoving volume of sources emitting CRs of type i. 
We now discuss the calculation of these terms and their scaling with redshift. 

The angular-averaged (differential) interaction rate, Ti (7^), appearing on the r.h.s. of Eq. (99) is 
defined as 



1 

Ti(z,Ei) = ^ J d cos 9 J duj{l- (3cos9)n 1 (z,uj)c tot 
-1 



, (100) 



ll3 (z,E u E 3 ) = Ti&Ed^iEitEj), (101) 

where n 7 (z,cj) is the energy distribution of background photons at redshift z and dNij/dEj is the 
angular- averaged distribution of particles j after interaction. The factor (1 — /3 cos 9) takes into account 
the relativistic Doppler shift of the photon density. 

In general, any transition i — >► i which can be approximated as 7^ (E, E') ~ 5(E—E f — AE)Ti(E) 
with AE/E <C 1 can be replaced in the Boltzmann equations (99) as 

-T(E)Y Z (E) + J dE' ju(E f , E) Yi(E f ) -+ d E (b z Y^ , (102) 

with bi = AETi ~ — E. The production of electron-positron pairs in the photon background with a 
small energy loss is usually approximated as a continuos energy loss (CEL) process. As we have seen 
in Sec. 2.3.1, it is also possible to approximate the energy loss in the hadronic cascade due to photopion 
production as a CEL 

^§(z,E) = bn(z,E) ~ ET p (z,E) - J dE 1 1 E' lpp (z,E,E') , (103) 

with &7r(0, E)/E given by (44) and (46). Diffractive processes at high energies with large final state 
multiplicities of neutrons and protons ultimately invalidate the CEL approximation. However, the relative 
error below 10 12 GeV is less than 15%, so we will use this approximation for a detailed numerical scan 



5 This is given by H 2 (z) — Hq [Q m (l + z) 3 + Q\], normalized to its value today of H ~ 100 h kms -1 Mpc -1 , in the 
usual "concordance model" dominated by a cosmological constant with Qa ~ 0.7 and a (cold) matter component, Q m ~ 
0.3 [15]. The time-dependence of the red-shift can be expressed via dz = —dt (1 + z)H. 
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Fig. 16: Left: The interaction and decay rates appearing in the Boltzmann equations for the CMB and CIB at 
z = 0. Right: Star formation rate (Eq. (Ill) from Ref. [158]) compared with our approximation (110) for the CIB 
photon density scaling with redshift. For comparison, the ex (1 + z) 3 scaling of the CMB photon density is also 
shown . From Ref. [157]. 



in the model space of proton spectra. For neutrons with energy less than 10 11 GeV the decay length is 
always smaller than the interaction length on the photon backgrounds. It is convenient to approximate 
the production of neutrons as ~ Y^ p + T£ n [157]; we have adopted this in all our calculations. In 
the left panel of Fig. (16) we show the quantities 6p air /E, 9^6 pair , T p and Ho for comparison. 

The redshift scaling of Eqs. (100) and (101) depend on how the photon backgrounds vary with 
redshift. The CMB spectral density (GeV -1 cm -3 ) scales adiabatically as: 

n 7 (z, u) = (1 + zf n 7 (0, u/(l + z)) , (104) 

following from Yy = dE(HEY 7 ) and Y 1 oc a 3 n 7 , where a is the cosmic scale factor. The scaling 
behavior of Eq. (104) translates into the following scaling of the quantities Yi and 7^ [21]: 

r t (z,E t ) = (i + z) 3 r,(o,(i + z)£;,), (105) 

ll3 {z,E u E 3 ) = (l + z)^ ij (0,(l + z)E i ,(l + z)E j ). (106) 

The scaling behaviour of b and its derivative c^fr do again depends on the photon background. For the 
CMB contribution, we have 

b^z.E,) = (l + z) 2 6,(0,(l + z)£;,), (107) 
d E b % {z,E % ) = (l + z ) 3 ^(0,(l + z)^). (108) 

We now discuss the red-shift scaling of the quantities bi, dsbi, jij and Ti for the case of the cosmic 
infrared background. The CIB spectrum has been studied and tabulated for redshifts up to z = 2 [159] 
(which we extrapolate slightly to UV energies as seen in Fig. 17). This is consistent with the constraints 
on the 7-ray opacity of the universe set by HESS [161], MAGIC [162] and Fermi-LAT [163] from non- 
observation of the expected cutoffs in the 7-ray spectra of extragalactic sources. The redshift dependence 
is given by: 

roc ^ 

ncmO*, (1 + z)E) = (1 + zf \ dz' ——Q C m(z' , (1 + z')E) , (109) 
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Fig. 17: The energy spectrum of the CMB and the CIB in the IR/optial and radio range at z = 0. The thin dashed 
line shows our extrapolation to UV energies. From Ref. [160]. 



where Qcib is the comoving luminosity density of the sources and we neglect absorption effects other 
than expansion. We assume that this follows the star formation rate: Qcm{z,E) oc HsfrX^) Qcib (0, £7). 
We can then infer the bolometric evolution to be: 

Ncmjz) _ 3 f™dz>Hs F K/(H(z>)(l + z>)) 

iVciB(O) U + f °°dz>HsFR/(H(z>)(l + z>))' U1U; 

where A^cib(^) is the number of infrared-optical photons per proper volume at redshift z. Following a 
recent compilation [158, 164], we adopt 



(l + zy A Z<1, 

7Vi(l + z)-°- 3 l<z<4, (HI) 
< 7Vi7V 4 (l + ^)" 3 - 5 z>4. 



with appropriate normalization factors, N\ = 2 3 7 and N4 = 5 3 2 (see the right panel of Fig. 16). For 
comparison, the CMB evolves as A^cmb(^) /^Vcmb(0) = (1 + z) 3 . To simplify the numerical evaluation 
we approximate the evolution with redshift as 

n clB (z,u) ~ ^- ^mM nc iB(0,a;/(l + *)) , (112) 

and this is shown in the right panel of Fig. 16. The redshift scaling of the quantities 7^ , b{ and 

for the CIB is then obtained from the corresponding scaling given for the CMB in Eqs. (105/106) and 

(107/108), by multiplying the r.h.s. by a factor N Qm {z) /Ncib(0)/(1 + zf . 

We have little direct knowledge of the cosmic radio background. An estimate made using the RAE 
satellite [165] is often used to calculate the cascading of UHE photons. A theoretical estimate has been 
made [166] of the intensity down to kHz frequencies, based on the observed luminosity function and 
radio spectra of normal galaxies and radio galaxies, although there are large uncertainties in the assumed 
evolution. The calculated values are about a factor of ^ 2 above the measurements and to ensure maximal 
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energy transfer in the cascade we will adopt this estimate and assume the same redshift scaling as the 
cosmic infrared/optical background. We summarize the adopted cosmic radiation backgrounds in Fig. 17. 

The emission rate of CR protons per comoving volume is assumed, as per usual practice, to follow 
a power-law: 

{/_ [E j i?min) E < E m [ n , 
1 £min<£<£max, (H3) 

f+(E/ £ max ) ^max < E . 

We will consider spectral indices 7 in the range 2-1-3. The functions f±(x) = x ±2 exp(l — x ±2 ) in 
Eq. (113) smoothly turn off the contribution below i? m in and above i? m ax- We set £ max = 10 12 GeV 
in the following and vary E m i n in the range 10 8 5 -7- 10 10 GeV, corresponding to a galactic-extragalactic 
crossover between the second knee and the ankle in the CR spectrum. 

The cosmic evolution of the spectral emission rate per comoving volume is parameterized as: 

Qp(z,E)=H(z)Qp(0,E). (114) 
For simplicity, we use the standard approximation 

n(z) = (l + z) n &(z m ^-z), (115) 

with z max = 2. Note that the dilution of the source density due to the Hubble expansion is taken care of 
since Q is the comoving density, i.e. for no evolution we would simply have H = 1. We consider cosmic 
evolution of UHECR sources with n in the range 2-1-6. 

We can express the system of partial integro-differential equations (99) in terms of a system of 
ordinary integro-differential equations [157], 

£ i = -HS i -b i (z,S i ), (116) 
Zi = [Pi(z,£i) - Ti(z,£i)] Z i + (l + z)Qf{z,£ i ) , (117) 

where we have defined f3i(z,E) = dEbi(z,E) and Zi(z,E) = (1 + z)Yi(z,£i(z,E)). The quantity 
£i(z, E) gives the energy that a particle of type i had at redshift z if we observe it today with energy E 
and take into account CEL. The effective source term in Eq. (1 17) is 

Qf{z,£ i {z,E)) = Q i + Y, J dEd E £ 3l3% {z,£ 3 ,E % )^-, (118) 

3 

where £j (z , E) and Zj (z, E) are subject to the boundary conditions £j (0, E) = E and Zj (^ m ax, E) = 0. 
The flux at z = can be expressed as 



ME) = ^Zi(0,E) 



^ />oo r pz 

— — dzexr> / dz' 
4vr J [Jo 



j d E b i (z',e i (z',E))-T i (z',ei(z',E)) 
(l + z')H(z>) 



H(z) 



For the numerical evaluation of the Boltzmann equations (99) it is convenient to solve [160] 

-H(z)(l + z)d z 3fi(z, E) = -T(z, £)m(z, E) + -^—^d E [b(z, £)%{z, E)\ + (1 + z)Qf(z, S) . 

(120) 

where we have defined £ — (l+z)E, and 3?i(z, E) = (l+z)Yi(z, $), subject to the boundary condition 
^K^max, E) = 0. The effective source term becomes Qf ff (z, <§) = Q% + J2j J dEj jji(z, $j,Ei) 3?j. 
In our calculation we use logarithmic energy bins with size Alog 10 E = 0.05 between 10 5 GeV and 
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10 5 GeV. For the numerical evolution of the differential Eq. (1 16) we use a simple Euler method with a 
step-size Az = 1CT 4 . The corresponding step-size in the propagation distance Ar = cAt is then always 
smaller than the proton interaction length. The flux at z = is given by Ji{E) = i^(0, E)/ (47r). 

Electromagnetic (EM) interactions of photons and leptons with the extra-galactic background light 
and magnetic field can happen on time-scales much shorter than their production rates. It is convenient 
to account for these contributions during the proton propagation as fast developing electro-magnetic cas- 
cades at a fixed redshift. We will use the efficient method of "matrix doubling" [167] for the calculation 
of the cascades. Since the cascade 7-ray flux is mainly in the GeV-TeV region and has an almost univer- 
sal shape here, it is numerically much more efficient to calculate the total energy density cj cas injected 
into the cascade and compare this value to the limit imposed by Fermi-LAT. The total energy density (eV 
cm -3 ) of EM radiation from proton propagation in the past is given as 

^ cas = J dEEncss^E) = JdtJ dE b ™^' z ^ np(z,E), (121) 

where n p (z, E) is the physical energy density at redshift z, defined via n p (z, E) = (1 + z) 3 Y p (z, E)\ 
details on the calculation are given in Appendix B. The continuous energy loss of protons into the 
cascade, denoted by 6 cas , is in the form of electron, positron and 7-ray production in BH (&bh) an d 
photopion (b n ) interactions. We derive the BH and photopion contribution to cj cas separately. For the 
photopion contribution we estimate 

b n (z,E) ~ J dE'E' [ lpe -(z,E,E') + lpe+ (z,E,E') + lpi (z,E,E')] , (122) 

where the angular-averaged distribution of particle j after the interaction are determined using the Monte 
Carlo package SOPHIA [168]. For the energy loss via BH pair production we use (37). Note, that since 
the photopion contribution in the cascade is dominated at the GZK cutoff, its contribution should increase 
relative to BH pair production with increasing crossover energy and, hence, also the associated neutrino 
fluxes after normalization to 7-ray and CR data. 

Next, we study the constraint set by the diffuse 7-ray background on all-proton models of ex- 
tragalactic CRs. We parametrize our ignorance of the crossover energy - which marks the transition 
between the galactic and extragalactic components - as a variable low energy cutoff in the proton in- 
jection rate. By fitting only to CR data above the crossover energy, taken to be between 10 8 5 GeV and 
10 10 GeV, we determine the statistically preferred values of the spectral index 7, cosmic source density 
evolution index n, and crossover energy E m [ n by a goodness-of-fit (GOF) test of the HiRes data, taking 
into account the energy resolution of about 25%. For each model we check that the total energy density 
of the EM cascade is below a critical value inferred from the recent measurement of the extragalactic 
7-ray background by the Fermi-LAT Collaboration [169]. 

Given the acceptance A{ (in units of area per unit time per unit solid angle) of the experiment for 
the energy bin i centered at E{ and with bin width A^, and the energy scale uncertainty of the experiment, 
(Te s the number of expected events in the bin is given by 

Ni{n^N^S) = Ai / J p Nqu (E)dE, (123) 

where 

J P N ,n,,(E) = ^n p (0,E) (124) 

is the proton flux arriving at the detector corresponding to a proton source luminosity as in Eq. (113), with 
the cosmic evolution of the source density given by Eqs. (114) and (115). The parameter S in Eq. (123) 
above is a fractional energy-scale shift that reflects the energy-scale uncertainty of the experiment, and 
Nq is the normalization of the proton source luminosity. 
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Fig. 18: Left Panel: Goodness-of-fit test of the HiRes data [28]. We show the 68% (pink), 95% (blue) and 99% 
(magenta) confidence levels of the injection index 7 and the cosmic evolution index n . The black lines indicate the 
allowed regions before the cascade (co> cas ) bound is imposed. Right Panel: The corresponding energy density in 
the EM cascade. From Ref. [160]. 



The probability distribution of events in the i-th bin is of the Poisson form with mean A^. Corre- 
spondingly the r-dimensional (r being the number of bins of the experiment with Ei > E^^) probability 
distribution for a set of non-negative integer numbers k = {fci, ...k r }, Pj^{n, 7, iVo, 5), is just the product 
of the individual Poisson distributions. According to this r-dimensional probability distribution, the ex- 
perimental result iV exp = {A^ xp , A^ xp } has a probability Pfi exp (n, 7, iVo, S) and, correspondingly, 
the experimental probability after marginalizing over the energy scale uncertainty and normalization is: 

(n,7,JV ,<5), (125) 

where the maximization is made within some prior for 6 and Nq. For the energy shift 6 we have used 
two forms for the prior, either a top-hat spanning the energy- scale uncertainty of the experiment, cfe s , or 
a gaussian prior of width <je s • 

For No we impose the prior arising from requiring consistency with the Fermi-LAT measure- 
ments [169] of the diffuse extra-galactic 7-ray background. In order to do so, we obtain the total energy 
density of EM radiation from the proton propagation using Eq. (121) and we require following Ref. [170]: 



,(N 0l n, 7) < 5.8 x 10" 7 eV/cm 3 . (126) 



The marginalization in Eq. (125) also determines iVQ est and S hest for the model, which are the 
values of the energy shift and normalization that yield the best description of the experimental CR data, 
subject to the constraint imposed by the Fermi-LAT measurement. 

Altogether the model is compatible with the experimental results at given GOF if 



J2 Pk(n, 7, < eSt , S best )0 \p n (n, 7, iV bes \ 5 best ) - P exp (n, 7) 



< GOF . (127) 



Technically, this is computed by generating a large number iV re p of replica experiments according to 
the probability distribution P%(n, 7, NQ est , S hest ) and imposing the fraction F of those which satisfy 

P^n, 7, JV , 6 hest ) > P exp (n, 7) to be F < GOF. 
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Fig. 19: The allowed proton flux (at the 99% confidence level) for increasing crossover energy E m { n . Each fit of 
the proton spectrum is marginalized with respect to the experimental energy uncertainty and we show the shifted 
predictions in comparison to the HiRes central values [28]. For comparison we also show the Auger data [29, 30] 
which has not been included in the fit. From Ref. [160]. 



With this method we determine the value of (n, 7) parameters that are compatible with the HiRes 
I and HiRes II experiments [28]. In the left panel of Fig. 18, we plot the regions with GOF 64%, 
95% and 99% for four values of the minimum (i.e. crossover) energy. In the right panel, we show 
the corresponding ranges of cj cas? best f° r the models as a function of the cosmic evolution index n. In 
order to display explicitly the impact of the constraint from the Fermi-LAT measurements of the diffuse 
extra-galactic 7-ray background (126), we show the corresponding GOF regions without imposing that 
constraint. In Table. 3 we list the parameters corresponding to the best-fit models and to the models 
with minimal and maximal contributions to uo^ and cj cas = uo^ + c^bh at the 99% C.L., together with 
the corresponding energy shifts which give best fits to the HiRes I and HiRes II data. We also show the 
parameters for the models with maximum uj^ and cj cas without imposition of the Fermi-LAT constraint. 

It is interesting to compare the allowed values of n with those in the cosmological evolution of 
UHECR candidate sources. For GRBs [171], 

^grb(^) = (1 + ^) L4 ^sfrW, (128) 

for AGNs [172,173] 



'HagnO) 



(l + z) 5 - u z < 1.7, 

N lm7 1.7<z<2.7, (129) 
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Fig. 20: Systematic effect of the experimental energy resolution on the fitted spectral index 7 and cosmological 
evolution parameter n. For illustration we show the dependence of the 95% C.L. bound for a crossover energy of 
10 18 eV. The blue contour corresponds to the region shown in Fig. 18 assuming an uncorrelated energy shift of 
25% in both data sets (HiRes I and II) [28], for a flat prior ("top-hat" distribution). The red dashed curve assumes 
correlated errors of the energy resolution in both data sets. The black dotted curve shows the result for uncorrelated 
errors with a Gaussian prior, and the dashed-dotted line shows uncorrelated errors with a flat prior, but with a lower 
uncertainty of 15%. From Ref. [160]. 



with TVi.r 2.7 5 and N 2 . 7 



10 043 , and for QSO [174] 

f(l + z) 3 
= < (1 + 1.9) 3 

[(l + 1.9) 3 exp{(2.7-z)/2.7}, 



z < 1.9, 

1.9 < z < 2.7, (130) 
z > 2.7. 



As an illustration of the agreement with the CR data, in Fig. 18 we show the range of proton 
fluxes corresponding to models with GOF 99% or better for increasing crossover energies E m [ n . As 
discussed above, each fit of the proton spectra is marginalized with respect to the experimental energy 
scale uncertainty and we show the shifted predictions with S hest in comparison to the HiRes data at 
central value. In the figure we also show the results from Auger [29, 30], though these have not been 
included in the analysis. 

These results are obtained assuming an energy scale uncertainty <je s = 5% with a "top-hat" prior 
for the corresponding energy shifts which are taken to be uncorrelated for HiRes I and HiRes II. In Fig. 20 
we explore the dependence of the results on these assumptions by using a different form for the prior, 
assuming the energy shifts to be correlated between the two experiments, or reducing the uncertainty to 
ce s = 15%. As seen in the figure, the main effect is associated with the reduction of the energy scale 
uncertainty which, as expected, results in a worsening of the GOF for models with larger n. This is 
directly related to the normalization constraint from Eq. (126). If one naively ignores the energy scale 
uncertainty, the constraint in Eq. (126) rules out models with n > 3 (the precise value depending on the 
assumed E m { n ). However, once the energy scale uncertainty is included, the constraint of Eq. (126) plays 
a weaker role on the determination of the GOF of the models. It does, however, imply a maximum value 
of A^Q est which, as we will see, impacts the corresponding ranges of neutrino fluxes. 
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Table 3: Cosmic ray source parameters which best fit the HiRes data [28], along with those which yield minimal 
and maximal contributions to uj^ (i.e. neutrino fluxes) and uj cas = uj^ -f co>bh (i.e. 7-ray fluxes), all at the 99% C.L. 







£ min = 10 8 - 5 GeV 


E m[n = 10 9 GeV 




mouei 


n 


7 


1 > a 

^cas 


^/best ^J/best 


n 7 


, , a 
^cas 


^/best 


^7/best 




fit with Fermi-LAT bound: 


best fit 


3.50 


2.49 


5.8 


0.005 0. 


3.20 2.52 


5.2 


0.050 


0.045 


min. uj cas 


4.50 


2.31 


4.4 


-0.235 -0.245 


2.25 2.47 


1.7 


-0.120 


-0.150 


max. u cas 


4.60 


2.36 


5.8 


-0.185 -0.175 


3.35 2.55 


5.8 


0.050 


0.060 


min. uj^ 


2.00 


2.67 


4.9 


0.215 0.235 


2.00 2.51 


1.8 


-0.070 


-0.095 


max. uj n 


4.80 


2.29 


5.8 


-0.220 -0.215 


5.10 2.29 


5.8 


-0.250 


-0.250 




fit without Fermi-LAT bound: 


max. uj cas 


4.45 


2.44 


15 


0.135 0.155 


5.25 2.36 


27 


0.205 


0.205 


max. 


4.80 


2.36 


14 


0.050 0.055 


5.30 2.35 


26 


0.190 


0.190 







E min = 10 9 - 5 GeV 


E min = 10 10 GeV 




model 


n 
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, , a 
^cas 


^/best ^//best 


n 7 


1 > a 

^cas 


^/best 


dllbest 




fit with Fermi-LAT bound: 


best fit 


4.05 


2.47 


5.8 


0.015 0.005 


4.60 2.50 


4.4 


-0.030 


-0.065 


min. uj cas 


2.00 


2.45 


1.4 


-0.050 -0.060 


2.00 2.88 


0.44 


-0.220 


-0.250 


max. u cas 


4.95 


2.37 


5.8 


-0.165 -0.160 


4.45 2.13 


5.8 


0.130 


0.090 


min. 


2.00 


2.63 


2.1 


0.075 0.070 


2.00 2.88 


0.44 


-0.220 


-0.250 


max. uj n 


5.35 


2.28 


5.8 


-0.240 -0.250 


4.40 2.10 


5.8 


0.145 


0.100 




fit without Fermi-LAT bound: 


max. uj cas 


6.00 


2.49 


30 


0.120 0.135 


6.00 2.14 


23 


0.250 


0.210 


max. 


6.00 


2.47 


29 


0.120 0.125 


6.00 2.10 


23 


0.250 


0.210 



a in units of 10" 7 eV/cm 3 



The corresponding range of 7-ray and cosmogenic neutrino fluxes (summed over flavor) is shown 
in Fig. 21 for models with minimal and maximal energy density at the 99% C.L. As expected, the 
maximum 7-ray fluxes are consistent with the Fermi-LAT data within the errors. For illustration, we also 
show as a dotted line the "naive" 7-ray limit E 2 J C8iS < ccj™| x /47rlog(TeV/GeV), corresponding to a 
7-ray flux in the GeV-TeV range which saturates the energy density (126). 

The cosmogenic neutrino and photon fluxes are also sensitive to the primary composition. For 
example, if the primaries are iron nuclei one would expect a considerably lower flux on both photons and 
neutrinos, rendering photon pile-up measurements less helpful in constraining cosmogenic fluxes. (See 
Appendix A for details). 

2.5.3 Upper limits on the cosmic neutrino flux 

High energy neutrino detection is one of the experimental challenges in particle astrophysics for the 
forthcoming years. It is widely believed that one of the most appropriate techniques for neutrino detection 
consists of detecting the Cherenkov light from muons or showers produced by the neutrino interactions 
in underground water or ice. For a recent review see e.g. [175]. This allows instrumentation of large 
enough volumes to compensate for both the low neutrino cross section and the low fluxes expected. 
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Fig. 21: Comparison of proton, neutrino and 7-ray fluxes for different crossover energies. We show the best- fit 
values (solid lines) as well as neutrino and 7-ray fluxes within the 99% C.L. with minimal and maximal energy 
density (dashed lines). The values of the corresponding model parameters can be found in Table. 3. The dotted line 
labeled "maximal cascade" indicates the approximate limit E 2 J cas < cu^f^/Air log(TeV/GeV), corresponding 
to a 7-ray flux in the GeV-TeV range saturating the energy density (126). The 7-ray fluxes are marginally consis- 
tent at the 99% C.L. with the highest energy measurements by Fermi-LAT. The contribution around 100 GeV is 
somewhat uncertain due to uncertainties in the cosmic infrared background. From Ref. [160]. 



There are several projects under way to build sufficiently large detectors to measure the expected signals 
from a variety of neutrino sources. The IceCube facility, deployed near the Amundsen-Scott station, is 
the largest neutrino telescope in the world [176]. It comprises a cubic-kilometer of ultra-clear ice about 
a mile below the South Pole surface, instrumented with long strings of sensitive photon detectors which 
record light produced when neutrinos interact in the ice. 

CR experiments, like Auger, provide a complementary technique for UHECv detection by search- 
ing for deeply-developing, large zenith angle (> 75°) showers [177]. At these large angles, hadron- 
induced showers traverse the equivalent of several atmospheres before reaching detectors at the ground. 
Beyond about 2 atmospheres, most of the electromagnetic component of a shower is extinguished and 
only very high energy muons survive. Consequently, a hadron-induced shower front is relatively flat and 
the shower particles arrive within a narrow time window (Fig. 22, top panel). In contrast, a neutrino 
shower exhibits characteristics similar to those of a vertical shower, which has a more curved front and a 
wider distribution in particle arrival times due to the large number of lower energy electrons and photons. 
Furthermore, the "early" part of the shower will tend to be dominated by the electromagnetic component, 
while "late" portion will be enriched with tightly bunched muons (Fig. 22, middle panel). Using these 
characteristic features, it is possible to distinguish neutrino induced events from background hadronic 
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showers. Moreover, because of full flavor mixing, tau neutrinos are expected to be as abundant as other 
species in the cosmic flux. Tau neutrinos can interact in the Earth's crust, producing r leptons which may 
decay above the Auger detectors [178-180] (Fig. 22, bottom panel). Details on how such events can be 
selected at the Auger Observatory are discussed in [181, 182]. 




Fig. 22: Schematic illustration of the properties of a hadron-induced shower (top), an z/-induced nearly horizontal 
shower (middle) and a ^-induced earth skimming shower (bottom). Note that only up-going showers resulting 
from t neutrino interactions in the Earth can be detected with any efficiency using the surface array. In contrast, 
all three neutrino species can be detected in down-going showers. Also note from the inset in the lower panel, that 
the incident r can experience several CC interactions and decays and thereby undergo a regeneration process. The 
Andes mountain range lies to the west of the observatory, and provides roughly an additional 20% target volume 
for v T interactions. 

So far, no neutrino candidates have been observed resulting in upper limits on the diffuse flux of 
neutrinos. For the case of up-going r neutrinos in the energy range 2 x 10 8 GeV < E v < 2 x 10 10 GeV, 
assuming a diffuse spectrum of the form E~ 2 , the current 90%CL bound, is [183] 

El ~ 4.7;^ x 10~ 8 GeV cm- 2 s" 1 sr" 1 . (131) 

Though Auger was not designed specifically as a neutrino detector, it is interesting to note that it exhibits 
good sensitivity in an energy regime complementary to those available to other dedicated instruments. 
For example, the 90% CL upper limit from IceCube, for neutrinos of all flavors in the energy range 
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2.0 x 10 6 GeV < E v < 6.3 x 10 9 GeV, is [184] 

E 2 ^i/aii _ 3 g x 1Q -8 Q e y cm -2 s -l gr -l _ ( j 32) 

We now derive model-independent bounds on the total neutrino flux. The event rate for quasi- 
horizontal deep showers is 

N = E / dEi Na °iN^x(Ei) £{E X ) , (133) 

i,x J 

where the sum is over all neutrino species i = v e , v e , v^, v^, v T , v T , and all final states X. Na = 
6.022 x 10 23 is Avogadro's number, and <fr l is the source flux of neutrino species i, a as usual denotes 
the cross section, and £ is the exposure measured in cm 3 w.e. sr time. The Pierre Auger Collaboration 
has searched for quasi-horizontal showers that are deeply-penetrating [183]. There are no events that 
unambiguously pass all the experimental cuts, with zero events expected from hadronic background. 
This implies an upper bound of 2.4 events at 90%CL from neutrino fluxes [185]. Note that if the number 
of events integrated over energy is bounded by 2.4, then it is certainly true bin by bin in energy. Thus, 
using Eq. (133) one obtains 

V / dE i N A <S> i (E i )a iN ^ x (E i )£(E i ) <2.4, (134) 

7^ Ja 

at 90% CL for some interval A. Here, the sum over X takes into account charge and neutral current 
processes. In a logarithmic interval A where a single power law approximation 

&(Ei) (JiN^xiEi) £(Ei) ~ Ef (135) 

is valid, a straightforward calculation shows that 

[ {E)e ' ^E^ a iN ^ x £ = (a iN ^ x £ E % **) ^ A , (136) 

J(E)e- A / 2 Ei 

where S = (a + l)A/2 and (A) denotes the quantity A evaluated at the center of the logarithmic 
interval [186]. The parameter a = 0.363 + f3 — 7, where 0.363 is the power law index of the SM 
neutrino cross section [187] and f3 and —7 are the power law indices (in the interval A) of the exposure 
and flux <3> z , respectively. Since sinhS/S > 1, a conservative bound may be obtained from Eqs. (134) 
and (136): 

N A ^2((TiN^x(Ei)) (S(Ei)) (Ei&) < 2.4/A . (137) 

i,X 

By taking A = 1 as a likely interval in which the single power law behavior is valid (this corresponds 
to one e-folding of energy), it is straightforward to obtain upper limits on the neutrino flux. The model- 
independent upper limits on the total neutrino flux, derived using an equivalent of 0.8 yr of full Auger 
exposure, are collected in Table 4 [183]. 



3 Phenomenology of extensive air showers 

3.1 Systematic uncertainties in air shower measurements from hadronic interaction models 

Uncertainties in hadronic interactions at ultrahigh energies constitute one of the most problematic sources 
of systematic error in the analysis of air showers. This section will explain the two principal schools of 
thought for extrapolating collider data to ultrahigh energies. 

Soft multiparticle production with small transverse momenta with respect to the collision axis is 
a dominant feature of most hadronic events at cm. energies 10 GeV < < 62 GeV (see e.g., [188, 
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Table 4: Model-independent upper limits on the neutrino flux at 90% CL. 



E v (GeV) (Ey gjg (cm" 2 sr" 1 s" 1 ) 

1 x 10 8 4.3 x 10- 14 

3 x 10 8 5.3 x 1(T 15 

1 x 10 9 1.2 x l(r 15 

3 x 10 9 4.7 x l(r 16 

1 x 10 10 2.2 x 10~ 16 

3 x 10 10 1.3 x 10- 16 

1 x 10 11 7.2 x 1CT 17 

3 x 10 11 4.3 x 1(T 17 



189]). Despite the fact that strict calculations based on ordinary QCD perturbation theory are not feasible, 
there are some phenomenological models that successfully take into account the main properties of the 
soft diffractive processes. These models, inspired by 1/N QCD expansion are also supplemented with 
generally accepted theoretical principles like duality, unitarity, Regge behavior, and parton structure. The 
interactions are no longer described by single particle exchange, but by highly complicated modes known 
as Reggeons. Up to about 62 GeV, the slow growth of the cross section with y/s is driven by a dominant 
contribution of a special Reggeon, the Pomeron. 

At higher energies, semihard (SH) interactions arising from the hard scattering of partons that 
carry only a very small fraction of the momenta of their parent hadrons can compete successfully with 
soft processes [190-197]. These semihard interactions lead to the "minijet" phenomenon, i.e. jets with 
transverse energy (Et = \p T \) much smaller than the total cm. energy. Such low-p T processes cannot 
be identified by jet finding algorithms, but (unlike soft processes) still they can be calculated using 
perturbative QCD. The cross section for SH interactions is described by 

"QCD(«,j£ in ) = £ / — [— d\t\ si/iOri, \t\) x 2 f j (x 2 , \t\) , (138) 

where x\ and X2 are the fractions of the momenta of the parent hadrons carried by the partons which col- 
lide, d&ij/d\i\ is the cross section for scattering of partons of types i and j according to elementary QCD 
diagrams, fi and fj are parton distribution functions (PDFs), s = x\ X28 and —i=s (1— cos #*)/2 = Q 2 
are the Mandelstam variables for this parton-parton process, and the sum is over all parton species. Here, 

p T = E^ sintf jet = ^ sinr , (139) 

and 

p„ =£] e a t b COS^jet, (140) 

where E 1 ^ is the energy of the jet in the lab frame, $j et the angle of the jet with respect to the beam 
direction in the lab frame, and #* is the angle of the jet with respect to the beam direction in the cm. 
frame of the elastic parton-parton collision. This implies that for small p 2 T « Q 2 . The integration 
limits satisfy 

Q 2 min < \t\ <s/2, (141) 

where Qmin = 1 — 2 GeV is the minimal momentum transfer. The measured minijet cross sections indi- 
cate that the onset of SH interactions has just occurred by CERN SPS energies (y/s > 200 GeV [198]. 

A first source of uncertainty in modeling cosmic ray interactions at ultrahigh energy is encoded in 
the extrapolation of the measured parton densities several orders of magnitude down to low x. Primary 
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protons that impact on the upper atmosphere with energy ~ 10 11 GeV yield partons with x = 2p* / ~ 
vn^j ~ 10 -7 , whereas current data on quark and gluon densities are only available for x > 10 -4 to 
within an experimental accuracy of 3% for Q 2 « 20 GeV 2 [199]. In Fig. 23 we show the region of the 
x — Q 2 plane probed by HI, ZEUS, and fixed target experiments. In addition, extrapolation of HERA 
data to UHECR interactions assumes universality of the PDFs. This assumption, based on the QCD 
factorization conjecture (the cross section of Eq. (138) can always be written in a form which factorizes 
the parton densities and the hard interaction processes irrespective of the order in perturbation theory 
and the particular hard process) holds in the limit Q 2 » Aqcd, where Aqcd ~ 200 MeV is the QCD 
renormalization scale. 

For large Q 2 and not too small x, the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equa- 
tions [201-204] 

9 (g^Q 2 )\ = a,(Q 2 ) fP qq P q9 \ - (q(x,Q 2 )\ (U2i 

d\nQ 2 \g(x,Q 2 )J 2n \P gq Pj ™ \g(x,Q 2 )J K } 

successfully predict the Q 2 dependence of the quark and gluon densities (q and respectively). Here, 
a s = g 2 /(4:7r), with g the strong coupling constant. The splitting functions indicate the probability 
of finding a daughter parton i in the parent parton j with a given fraction of parton j momentum. This 
probability will depend on the number of splittings allowed in the approximation. In the double-leading- 
logarithmic approximation, that is lim^^o ln(l/x) and limg2^ 00 1ii((5 2 /Aqcd), the DGLAP equations 
predict a steeply rising gluon density, xg ~ x -cu , which dominates the quark density at low x. This 
prediction is in agreement with the experimental results from HERA shown in Fig. 24 [206,207]. HERA 
data are found to be consistent with a power law, xg(x, Q 2 ) ~ x~ Ali , with an exponent Ah between 0.3 
and 0.4 [208]. 
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Fig. 24: Gluon momentum distributions xg(x, Q 2 ) in the proton as measured by the ZEUS and HI experiments at 
various Q 2 . From Ref. [205]. 



The high energy minijet cross section is then determined by the dominant gluon distribution 

^qcd^pT) « f /2 d\t\ ^| x ig ( Xl ,\i\) x 2 g(x 2 ,\t\) , (143) 

J xi J x 2 Jq2. d\t\ 

where the integration limits satisfy 



x\x 2 s > 2\t\ > 2Q 



2 

min ' 



(144) 



Furthermore, because da/d\i\ is peaked at the low end of the \i\ integration (see e.g [209]), the high 
energy behavior of ctqcd is controlled (via the lower limits of the x\, x 2 integrations) by the small-x 
behavior of the gluons [210] 



( \ f dx ^ -A H f dx 2 -A H A h i 

J2Q 2 . Is x l J2Q 2 . Ixi8 x< 2 



s 



(145) 



53 



This estimate is, of course, too simplistic. At sufficiently small x, gluon shadowing corrections suppress 
the singular x~ Ah behavior of xg and hence suppress the power growth of ctqcd with increasing s. 

Although we have shown that the onset of semihard processes is an unambiguous prediction of 
QCD, in practice it is difficult to isolate these contributions from the soft interactions. Experimental 
evidence indicates that SH interactions can essentially be neglected up to and throughout the CERN ISR 
energy regime, y/s < 62 GeV. Therefore, measurements made in this energy region can be used to 
model the soft interactions. A reasonable approach introduced in [21 1] is to base the extrapolation of the 
soft interactions on the assumption of geometrical scaling [212], which is observed to be true throughout 
the ISR energy range [213,214]. To this end, we introduce the standard partial- wave amplitude in impact- 
parameter space f(s,b), which is the Fourier transform of the elastic pp (or pp) scattering amplitude. (We 
neglect any difference between pp and pp for > 200 GeV.) Geometrical scaling (GS) corresponds 
to the assumption that /, which a priori is a function of two dimensional variables b and s, depends only 
upon one dimensional variable (3 = b/R(s), where R is the energy dependent radius of the proton, i.e. 

f(s,b) = f GS (/3 = b/R(s)). (146) 

Physically, this means that the opaqueness of the proton remains constant with rising energy and that the 
increase of the total cross section, <T to t> in the ISR energy range reflects a steady growth of the radius 
R(s). An immediate obvious consequence of GS is that the partial wave at b = should be independent 
of energy 

/M = 0) = /gs(/? = 0). (147) 

Another consequence is that the ratio of elastic scattering to total cross section, a e \(s)/a to t(s), should 
be energy independent. This follows from 

otot = 8tt J lmf(s,b)bdb 

= 8tvR\s) Jlmf GS (P)Pdp, (148) 

GS 

and 

<Tel = 8^ J \f(s,b)\ 2 bdb 

= 8ttR 2 (s) J |/ GS (/3)| 2 /3d/3. (149) 

GS 

To determine the gross features at high energies we can assume that the elastic amplitude has a simple 
form 

F(s,t) = ia tot (s) e Bt /\ (150) 

with B the slope parameter that measures the size of the proton [215]. This is a reasonable assumption: 
the amplitude is predominantly imaginary, and the exponential behavior observed for \t\ < 0.5 GeV 2 
gives the bulk of the elastic cross section. Now, the Fourier transform f(s,b) of the elastic amplitude 
t) given by (150) has a Gaussian form in impact parameter space 

/M) = ^, (151) 

and it follows that 

Im/( S ,6 = 0) = ^| = 2 ^. (152) 

87T.O cr tot 
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Equation (152) offers a very clear way to see the breakdown of GS and to identify semihard interactions 
from the growth of the central partial wave. 

In general unitarity requires Im/(s, 6) < \, which in turm implies cr e \/(Jtot < \ [215]. This seems 
to indicate that the Gaussian form (151) may not longer be applicable at ultrahigh energies, but rather it 
is expected that the proton will approximate a "black disk" of radius bo, i.e. /(s, b) = | for < b < bo 
and zero for b>bo. Then a e \ ~ ^cr to t — ^b\. 

In order to satisfy the unitarity constraints, it is convenient to introduce the eikonal x defined by 

f( s ,b)= t -{l-exp[i x (s,b)]} , (153) 

where Imx > 0. If we neglect for the moment the shadowing corrections to the PDFs and take xg oc 
X -A H j n ^ sma n_ x ii m it W e obtain, as explained above, power growth of the cross section for SH 
interactions, oqcd ~ s Au and Imx(s,6 = 0) » 1 as s — >► oc. Indeed we expect Imx ^ 1 ( an d 
unitarity to be saturated) for a range of b about 6 = 0. Then we have 

roc 

<rtot = 4^ / bdbQ(b -b) 
Jo 

rbo(s) 

~ 4tt / bdb = 2irbl, (154) 
Jo 

with x — Xsh an( i where bo(s) is such that 

Mo(*))^l. (155) 



Hereafter, we ignore the small real part of the scattering amplitude, which is good approximation 
at high energies. The unitarized elastic, inelastic, and total cross sections (considering now a real eikonal 
function) are given by [216-219] 

cr e i = 27r J dbb {1 -exp[-x soft (5,6) -x SH ( 5 ^)]} 2 > ( 156 ) 
= 2n J dbb{l-exp[-2x soit (s,b)-2x Sii (s,b)}} , (157) 
cr t ot = 47r y dbb {1 -exp[-x soft (s,6) -X SH ( 5 ^)]} ? 



Cinel 



(158) 

where the scattering is compounded as a sum of QCD ladders via SH and soft processes through the 
eikonals x SH and x soft - 

Now, if the eikonal function, x(s, b) = x soft (<§, 6) + Xsh ( s ^) = indicates the mean number 
of partonic interaction pairs at impact parameter 6, the probability p n for having n independent partonic 
interactions using Poisson statistics reads, p n = (X n /n\) e~ x . Therefore, the factor 1 — e~ 2x = Yl^Li Pn 
in Eq. (157) can be interpreted semiclassically as the probability that at least 1 of the 2 protons is broken 
up in a collision at impact parameter b. With this in mind, the inelastic cross section is simply the 
integral over all collision impact parameters of the probability of having at least 1 interaction, yielding 
a mean minijet multiplicity of (n m i n ij et ) ~ CQCD/^inei [220]. The leading contenders to approximate 
the (unknown) cross sections at cosmic ray energies, SIBYLL [221] and QGSJET [222], share the eikonal 
approximation but differ in their ansdtse for the eikonals. In both cases, the core of dominant scattering 
at very high energies is the SH cross section given in Eq. (138), 

X S H = l*QCD(s,p™ in )A(s,b), (159) 
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where the normalized profile function, 2tt dbbA{s 1 b) = 1, indicates the distribution of partons in 
the plane transverse to the collision axis. 

In the QGS JET-like models, the core of the SH eikonal is dressed with a soft-pomeron pre-evolution 
factor. This amounts to taking a parton distribution which is Gaussian in the transverse coordinate dis- 
tance 6, 

e -b 2 /R 2 (s) 

A(s,b)= , (160) 

7rR Z {S) 

with R being a parameter. For a QCD cross section dependence, oqcd ~ s Ah > one gets for a Gaussian 
profile 

b 2 (s) -i? 2 A H Ins (161) 

and at high energy 

rbo(s) 

a ine i = 2tt / d& 6 - Trir A H In 5 . (162) 
Jo 



If the effective radius R (which controls parton shadowing) is energy-independent, the cross section 
increases only logarithmically with rising energy. However, the parameter R itself depends on the colli- 
sion energy through a convolution with the parton momentum fractions, R 2 (s) ~ Rq + 4 a f eS In 2 5, with 
a' ef £ « 0.11 GeV -2 [145]. Thus, the QGSJET cross section exhibits a faster than In s rise, 

^inei ~ 4tt ag ff A H In 2 5 . (163) 

In SlBYLL-like models, the transverse density distribution is taken as the Fourier transform of the 
proton electric form factor, resulting in an energy-independent exponential (rather than Gaussian) fall-off 
of the parton density profile for large b, 

yb7T 

where K 3 denotes the modified Bessel function of the third kind and /j 2 « 0.71GeV 2 [221]. Thus, (159) 
becomes 

XsH ~e-^ S AH , (165) 

and (155) is satisfied when 

A H 

b (s) = — In s. (166) 

Therefore, for SlBYLL-like models, the growth of the inelastic cross section also saturates the In 2 s 
Froissart bound [223], but with a multiplicative constant which is larger than the one in QGS JET-like 
models 

A 2 

C^inel ~ 7TC —f In 2 S , (167) 

where the coefficient c « 2.5 is found numerically [145]. 

The main characteristics of the pp cascade spectrum resulting from these choices are readily pre- 
dictable: the harder form of the SIBYLL form factor allows a greater retention of energy by the leading 
particle, and hence less available for the ensuing shower. Consequently, on average SlBYLL-like models 
predict a smaller multiplicity than QGS JET-like models [42]. 

There are three event generators, sibyll [221], qgsjet [222], and dpmjet [224] which are tai- 
lored specifically for simulation of hadronic interactions up to the highest cosmic ray energies. The latest 
versions of these packages are sibyll 2.1 [225], qgsjet 11-03 [226], and dpmjet III [227]; respec- 
tively. In QGSJET, both the soft and hard processes are formulated in terms of Pomeron exchanges. To 
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describe the minijets, the soft Pomeron mutates into a "semihard Pomeron", an ordinary soft Pomeron 
with the middle piece replaced by a QCD parton ladder, as sketched in the previous paragraph. This 
is generally referred to as the "quasi-eikonal" model. In contrast, SIBYLL and DPMJET follow a "two 
channel" eikonal model, where the soft and the semi-hard regimes are demarcated by a sharp cut in the 
transverse momentum: SIBYLL 2.1 uses a cutoff parametrization inspired in the double leading logarith- 
mic approximation of the DGLAP equations, 

p™ in (^) = p ° T + 0.065 GeV exp[0.9 Vh^} , (168) 

whereas dpmjet III uses an ad hoc parametrization for the transverse momentum cutoff 

^ in (^) = P° T + 0.12 GeV [log 10 (^/50GeV)] 3 , (169) 

where p° T = 2.5 GeV [208]. 

The transition process from asymptotically free partons to colour-neutral hadrons is described in 
all codes by string fragmentation models [228]. Different choices of fragmentation functions can lead 
to some differences in the hadron multiplicities. However, the main difference in the predictions of 
QGSJET-like and SlBYLL-like models arises from different assumptions in extrapolation of the parton 
distribution function to low energy. 

The proton-air collisions of interest for air shower development cause additional headaches for 
event generators. Both SIBYLL and QGSJET adopt the Glauber formalism [216], which is equivalent to 
the eikonal approximation in nucleon-nucleon scattering, except that the nucleon density functions of the 
target nucleus are folded with that of the nucleon. The inelastic and production cross sections read: 

<ef ir ~ dbb i 1 - ex P K>t AT N (b)]} , (170) 

<"od r ~^Jdbb{l- exp [a inel AT N (b)}} , (171) 

where T/v(6) is the transverse distribution function of a nucleon inside a nucleus. Here, cr ine i and <7 to t 
are given by Eqs. (157) and (158), respectively. The p-air inelastic cross section is the sum of the "quasi- 
elastic" cross section, which corresponds to cases where the target nucleus breaks up without production 
of any new particles, and the production cross section, in which at least one new particle is generated. 
Clearly the development of EASs is mainly sensitive to the production cross section. Overall, the geo- 
metrically large size of nitrogen and oxygen nuclei dominates the inclusive proton-target cross section, 
and as a result the disagreement from model-dependent extrapolation is not more than about 15%. More 
complex nucleus-nucleus interactions are discussed in [229]. 



EXERCISE 3.1 Consider a typical air nuclei of average (A) = 14.5 and calculate the proton-air 
cross section using the approximated expressions (163) and (167) together with the z-integrated Woods- 
Saxon profile [230] 



T N (b) = | J°° dz {l + exp [(Vb 2 + z 2 - R N )/a\ } 1 , 



where 



47T 3 

z y N 



(172) 



(173) 



a — 0.5 fm and Rn — 1.1A 1 / 3 fm [231]. Compare the results with those shown in Fig. 25. 
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Table 5: Different k- values used in cosmic ray experiments. 



Experiment 


k 


Fly's Eye 


1.6 


Akeno 


1.5 


Yakutsk-99 


1.4 


EAS-TOP 


1.15 



Adding a greater challenge to the determination of the proton air cross section at ultrahigh en- 
ergies is the lack of direct measuremnts in a controlled laboratory environment. The measured shower 
attenuation length, A m , is not only sensitive to the interaction length of the protons in the atmosphere, 

Ap-air, With 

Am = fcA p _ air = fc^^, (174) 

°prod 

but also depends on the rate at which the energy of the primary proton is dissipated into EM shower 
energy observed in the experiment. Here, A m and A p _ a i r are in gem -2 , the proton mass m p is in g, and 
the inelastic production cross section cTp~^ ir is in mb. The value of k depends critically on the inclusive 
particle production cross section and its energy dependence in nucleon and meson interactions on the 
light nuclear target of the atmosphere. The measured depth X max at which a shower reaches maximum 
development in the atmosphere has been the basis of cross section measurements from experiments prior 
to HiRes and Auger. However, A max is a combined measure of the depth of the first interaction (which 
is determined by the inelastic cross section) and of the subsequent shower development (which has to be 
corrected for). The model dependent rate of shower development and its fluctuations are the origin of 
the deviation of k from unity in Eq. (174). As can be seen in Table 5, there is a large range of k values 
(from 1.6 for a very old model where the inclusive cross section exhibited Feynman scaling, to 1.15 for 
modern models with large scaling violations) that make the published values of cr p _ a i r unreliable. 

Recently, the HiRes Collaboration developed a quasi-model-free method of measuring o^^ 1 di- 
rectly [232]. This is accomplished by folding a randomly generated exponential distribution of first 
interaction points into the shower development program, and therefore fitting the entire distribution and 
not just the trailing edge. Interestingly, the measured k = 1.21^09 by tne HiRes group is in agreement 
with the one obtained by tuning the data to the theory [233, 234]. 

A compilation of published proton-air cross section measurements is shown in Fig. 25. In the 
left panel we show the data without any modification. In the right panel, the published values of o^~^ T 
for Fly's Eye [242], Akeno [238], Yakutsk-99 [240], and EAS-TOP [236] collaborations have been 
renormalized using the common value of k = 1.264 ± 0.033(stat) ± 0.013(syst) [234]. We have 
parametrized the rise of the cross section using a functional form that saturates the Froissart bound, 

a P; Q f = A- £ln(£/GeV) + Cln 2 (£/GeV) mb . (175) 

The curve with a fast rise, hereafter referred to as case-i, corresponds to A = 280, 5 = 5.7, and C = 0.9. 
The slow rise of case-iz has the following parameters: A = 290, B = 6.2, and C = 0.64. The behavior 
of the cross section in case-i roughly matches the one implemented in SIBYLL-2.1 [225]. 

Some guidance towards understanding hadronic processes in the forward direction may come 
directly from measurements of hadrons in airshowers [243]. However, the most useful experimental 
input in the forseeable future will likely come from the LHC. Processes with low momentum transfer 
tend to populate the region at very small angles $ with respect to the beam direction. The distribution 
of pseudorapidity, 77 = — lntan($/2), and the energy flow distribution are shown in Fig. 26. While the 
particle multiplicity is greatest in the low \rj\ region, it is clearly seen that the energy flow is peaked at 
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Fig. 25: Compilation of proton-air production cross section from cosmic ray measurements (ARGO-YBJ [235], 
EAS-TOP [236], HiRes [232], KASCADE prototype [237], Akeno [238], Yakutsk-90 [239], Yakutsk-99 [240], 
Tien-Shan EAS complex [241] and Fly's Eye [242]). The data are compared to the parametrizations discussed in 
the text; case-i corresponds to the dashed line and case-n to the dot-dashed line. 



small production angles (large \rj\). The two general-purpose experiments, ATLAS and CMS, cover up to 
\rj\ < 5. TOTEM's coverage in the very forward rapidity range, 3.1 < 77 < 6.5, is committed to measure 
the pp total and elastic scattering cross sections, see Appendix C. The fragmentation region that plays 
a crucial role in the development of EASs tends to populate the region corresponding to pseudorapidity 
range 6 < \rj\ < 10. 

The first batch of data collected by the ATLAS experiment has been used to study the rise of the 
pp cross section [245]. In Fig. 27 we show a comprehensive investigation into the physics underpin- 
ning the rise of cr to t and cr- ine \. 6 The figure shows several models which employ different assumptions 
about the evolution of the PDFs at low-x (GRV [249], GRV94 [250], GRV98 [251], MRST [252], and 
CTEQ [253]). Model I is an eikonal minijet model incorporating effects induced by the transverse 
momentum distribution of soft gluons [254]. In the spirit of the Bloch-Nordsieck study in electrody- 
namics [255], the parton distribution in b- space is determined as the Fourier transform of the resummed 
distribution of soft gluons (down to zero gluon momenta) emitted from the initial state during the colli- 
sion. The very large b limit of the profile function 

A(b,s) ~e-( A6 ) 2/ \ (176) 

yields an impact parameter distribution falling at its faster like a Gaussian ((3 = 1) and at its slowest like 
an exponential ((3 = 1/2). The scale A oc Aqcd includes a mild energy dependence as well as a residual 
dependence upon the parameter (3. This behavior in impact parameter space together with the high energy 
behavior of the minijet cross section leads to an asymptotic rise of the total cross section [256] 

* tot ~f 2 A^ InV^, (177) 

consistent with the Froissart bound. In Model II the effective density overlaping parton distributions in 
the colliding protons is a SlBYLL-like exponentially decreasing hyperbolic Bessel function, determined 

6 The CMS measurements of the pp cross section at y/s = 7 TeV (not shown in the figure) are summarized in [248]. 



59 



I§" 8 
O 

? 7 

6 
5 
4 
3 
2 
1 



W 1.4 

1 .2 



0.8 
0.6 
0.4 
0.2 




I i i i i I 



< 



LHC, inelastic collisions 

TOTEM, CMS 



Tracking and Calorimetry 



ATLAS, CMS 

< 7-7—, * ; > 



> 



Calorimetry 



ATLAS, CMS 
Tracking 



iii 



j— i i i i i — i — \ i i i i i i i i i i i i i i i i i i 



-2.5 -5 -2.5 2.5 5 2.5 10 



Fig. 26: Pseudorapidity distributions of charged particles (upper panel) and of the energy flow (lower panel) for 
pp collisions at LHC. From Ref. [244]. 



from the proton electric form factor [257]. The analyticity-constrained analytic amplitude model of 
Block and Halzen [258, 259] exploits high quality low energy data to determine an asymptotic growth 
of the cross section, <7 to t ~ In 2 s, which saturates the Froissart bound. The Aspen model is a revised 
version of the eikonal model of Block, Gregores, Halzen, and Pancheri [260], which now incorporates 
analyticity constraints [261]. 

In a complementary study, a comparison of LHC data (on inclusive particle production at ener- 
gies = 0.9, 2.36, and 7 TeV) with predictions of a variety of hadronic interaction models has been 
elaborated in [262]. The results show that while reasonable qualitative agreement has been achieved for 
some event generators, none of them reproduces the evolution of all the observables with particularly 
compeling precision. 

In summary, high energy hadronic interaction models are still being refined and therefore the 
disparity between them can vary even from version to version. At the end of the day, however, the relevant 
parameters boil down to two: the mean free path, AcR-air = (^p r od n atm) _1 , and the inelasticity, 
^cR-air = 1 — ^lead/^proj. where n atm is the number density of atmospheric target nucleons, £i e ad is the 
energy of the most energetic hadron with a long lifetime, and E pro j is the energy of the projectile particle. 
The first parameter characterizes the frequency of interactions, whereas the second one quantifies the 
energy lost per collision. Overall, SIBYLL has a shorter mean free path and a smaller inelasticity than 
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Fig. 27: Compilation of total and inelastic pp and pp cross section measurements. The data are compared to model 
predictions. The dashed lines indicate the predictions of Model II [246]. From low to high energy the same eikonal 
function has been used to compute both the total and the inelastic cross sections. The predictions are calculated 
using GRV pfd's, with p™ in = 1.15 GeV. The green band is obtained with Model I varying by no more 
than a few percent and using different partonic densities [246]. The blue band corresponds to Model II using 
pf in = 1.1 GeV (with MRST) in the upper edge and p™ in = 1.15 GeV (with GRV94) in the lower edge [246]. 
Predictions for a tot from the analytic amplitude model and for cr inel from the Aspen model are indicated by blue 
stars [247]. This figure is courtesy of Giulia Pancheri. 

QGSJET. Since a shorter mean free path tends to compensate for a smaller inelasticity, the two codes 
generate similar predictions for an air shower which has lived through several generations. Both models 
predict the same multiplicity below about 10 7 GeV, but the predictions diverge above that energy. Such 
a divergence readily increases with rising energy. While QGSJET predicts a power law-like increase of 
the number of secondaries up to the highest energy, SIBYLL multiplicity exhibits a logarithmic growth. 
As it is extremely difficult to observe the first interactions experimentally, it is not straightforward to 
determine which model is closer to reality. 

3.2 Electromagnetic processes 

The evolution of an extensive air shower is dominated by electromagnetic processes. The interaction of 
a baryonic cosmic ray with an air nucleus high in the atmosphere leads to a cascade of secondary mesons 
and nucleons. The first few generations of charged pions interact again, producing a hadronic core, which 
continues to feed the electromagnetic and muonic components of the showers. Up to about 50 km above 
sea level, the density of atmospheric target nucleons is n atm ~ 10 20 cm -3 , and so even for relatively 
low energies, say E n ± « 1 TeV, the probability of decay before interaction falls below 10%. Ultimately, 
the electromagnetic cascade dissipates around 90% of the primary particle's energy, and hence the total 
number of electromagnetic particles is very nearly proportional to the shower energy [263]. 

Roughly speaking, at 10 11 GeV, baryons and charged pions have interaction lengths of the order of 
40 g/cm 2 , increasing to about 60 g/cm 2 at 10 7 GeV. Additionally, below 10 10 GeV, photons, electrons, 
and positrons have mean interaction lengths of 37 g/cm 2 . Altogether, the atmosphere acts as a natural 
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colorimeter with variable density, providing a vertical thickness of 26 radiation lengths and about 15 
interaction lengths. Amusingly, this is not too different from the number of radiation and interaction 
lengths at the LHC detectors. For example, the CMS electromagnetic calorimeter is > 25 radiation 
lengths deep, and the hadron calorimeter constitutes 1 1 interaction lengths. 

By the time a vertically incident 10 11 GeV proton shower reaches the ground, there are about 10 11 
secondaries with energy above 90 keV in the the annular region extending 8 m to 8 km from the shower 
core. Of these, 99% are photons, electrons, and positrons, with a typical ratio of 7 to e + e~ of 9 to 1. 
Their mean energy is around 10 MeV and they transport 85% of the total energy at ground level. Of 
course, photon-induced showers are even more dominated by the electromagnetic channel, as the only 
significant muon generation mechanism in this case is the decay of charged pions and kaons produced in 
7-air interactions [264]. 

It is worth mentioning that these figures dramatically change for the case of very inclined showers. 
For a primary zenith angle, 9 > 70° , the electromagnetic component becomes attenuated exponentially 
with atmospheric depth, being almost completely absorbed at ground level. Note that the vertical atmo- 
sphere is « 1, 000 g/cm 2 , and is about 36 times deeper for completely horizontal showers. As a result, 
most of the energy at ground level from an inclined shower is carried by muons. 

In contrast to hadronic collisions, the electromagnetic interactions of shower particles can be cal- 
culated very accurately from quantum electrodynamics. Electromagnetic interactions are thus not a 
major source of systematic errors in shower simulations. The first comprehensive treatment of electro- 
magnetic showers was elaborated by Rossi and Greisen [265]. This treatment was recently cast in a more 
pedagogical form by Gaisser [63], which we summarize in the subsequent paragraphs. 

The generation of the electromagnetic component is driven by electron bremsstrahlung and pair 
production [266]. Eventually the average energy per particle drops below a critical energy, eo, at which 
point ionization takes over from bremsstrahlung and pair production as the dominant energy loss mech- 
anism. The e ± energy loss rate due to bremsstrahlung radiation is nearly proportional to their energy, 
whereas the ionization loss rate varies only logarithmically with the e ± energy. Though several differ- 
ent definitions of the critical energy appear in the literature [267], throughout these lectures we take the 
critical energy to be that at which the ionization loss per radiation lenght is equal to the electron energy, 
yielding e = 710 MeV/(Z eff + 0.92) ~ 86 MeV [268]. 7 The changeover from radiation losses to ion- 
ization losses depopulates the shower. One can thus categorize the shower development in three phases: 
the growth phase, in which all the particles have energy > eo; the shower maximum, X max ; and the 
shower tail, where the particles only lose energy, get absorbed or decay. 

The relevant quantities participating in the development of the electromagnetic cascade are the 
probability for an electron of energy E to radiate a photon of energy k = y hrem E and the probability 
for a photon to produce a pair e + e~ in which one of the particles (hereafter e~) has energy E = y pair k. 
These probabilities are determined by the properties of the air and the cross sections of the two processes. 

In the energy range of interest, the impact parameter of the electron or photon is larger than an 
atomic radius, so the nuclear field is screened by its electron cloud. In the case of complete screen- 
ing, where the momentum transfer is small, the cross section for bremsstrahlung can be approximated 



where A e $ is the effective mass number of the air, X EM is a constant, and is Avogadro's number. In 
the infrared limit (i.e. y hrem <C 1) this approximation is inaccurate at the level of about 2.5%, which is 
small compared to typical experimental errors associated with cosmic air shower detectors. Of course, 

7 For altitudes up to 90 km above sea level, the air is a mixture of 78.09% of N2, 20.95% of O2, and 0.96% of other 
gases [269]. Such a mixture is generally modeled as an homogeneous substance with atomic charge and mass numbers Z e ff = 
7.3 and A e s = 14.6, respectively. 



by [270] 




(178) 
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the approximation fails as y hvem — » 1, when nuclear screening becomes incomplete, and as y hrem — >> 0, at 
which point the LPM and dielectric suppression effects become important, as we discuss below. 

Using similar approximations, the cross section for pair production can be written as [270] 

<fo 7 -> e + e - _ Aeff ( 4 4 2 \ 

d# ~ x em tv a v 3?/pair 3 ^ pair /' 

The similarities between this expression and Eq. (178) are to be expected, as the Feynman diagrams for 
pair production and bremsstrahlung are variants of one another. 

The probability for an electron to radiate a photon with energy in the range (fc, k+dk) in traversing 
dt = dX/X BM of atmosphere is 

da e ^ X EM N A A+ _( ni 41- 2/ brem 



whereas the corresponding probability density for a photon producing a pair, with electron energy in the 
range (E, E + dE), is 

da 7 ^ e+e - X EM 7V A / 4 4 2 



« 1 - -y . + . ) dy . dt . (181) 

dE A eff V 3 3 pai V 

The total probability for pair production per unit of X EM follows from integration of Eq. (181), 

7 d# A eff y V 3 pair 3^ air y ^ pair 9 v ; 

As can be seen from Eq. (180), the total probability for bremsstrahlung radiation is logarithmically 
divergent. However, this infrared divergence is eliminated by the interference of bremsstrahlung ampli- 
tudes from multiple scattering centers. This collective effect of the electric potential of several atoms 
is known as the LPM effect [32, 33]. Of course, the LPM suppression of the cross section results in an 
effective increase of the mean free path of electrons and photons. This effectively retards the develop- 
ment of the electromagnetic component of the shower. It is natural to introduce an energy scale, ^lpm ? 
at which the inelasticity is low enough that the LPM effect becomes significant [271]. Below £?lpm, the 
energy loss rate due to bremsstrahlung is roughly 



dE 
dX 



f 2/ brem E L m + \ 1 J hre A dy hrem =-^~. (183) 

A EM ^0 V 6 Vbrem J A EM 



With this in mind, we now identify the constant X EM « 36.7 g cm - ^ with the radiation length in air 
defined as the mean distance over which a high-energy electron loses 1 je of its energy, or equivalently 
7/9 of the mean free path for pair production by a high-energy photon [267]. 

The most evident signatures of the LPM effect on shower development are a shift in the position 
of the shower maximum X max and larger fluctuations in the shower development. When considering 
the LPM effect in the development of air showers produced by UHECRs, one has to keep in mind that 
the suppression in the cross sections is strongly dependent on the atmospheric depth. 8 Since the upper 
atmosphere is very thin, the LPM effect becomes noticeable only for photons and electrons with energies 
above £^lpm ~ 10 10 GeV. For baryonic primaries the LPM effect does not become important until the 
primary energy exceeds 10 12 GeV. This is because the electromagnetic shower does not commence until 
after a significant fraction of the primary energy has been dissipated through hadronic interactions. To 
give a visual impression of how the LPM effect slows down the initial growth of high energy photon- 
induced showers, we show the average longitudinal shower development of 10 10 GeV proton and 7-ray 
showers (generated using AIRES 2.6.0 [273]) with and without the LPM effect in Fig. 28. 



8 The same occurs for dielectric suppression, although the influence is not as important as for the LPM effect [272]. 
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Fig. 28: Average longitudinal shower developments of 10 11 GeV proton (dashed-dotted line) and 7-rays with and 
without the LPM effect (solid and dotted lines, respectively). The primary zenith angle was set to 6 = 60°. The 
shadow area represents the intrinsic fluctuations of the showers. Larger fluctuations can be observed for 7-ray 
showers with the LPM effect, as expected. From Ref . [200] . 

At energies at which the LPM effect is important {viz. E > ElpmX 7-ray showers will have al- 
ready commenced in the geomagnetic field at almost all latitudes: primary photons with E > 10 10 GeV 
convert into e + e~ pairs, which in turn emit synchrotron photons. This reduces the energies of the pri- 
maries that reach the atmosphere, and thereby compensates for the tendency of the LPM effect to retard 
the shower development. The relevant parameter to determine both conversion probability and syn- 
chrotron emission is E x B±, where E is the 7-ray energy and B± the transverse magnetic field. This 
leads to a large directional and geographical dependence of shower observables. Thus, each experiment 
has its own preferred direction for identifying primary gamma rays. For instance, Fig. 29 shows a map 
of the photon conversion probability in the geomagnetic field for all incident directions evaluated at the 
location of the HiRes experiment (\B\ = 0.53 G, 1 = 25°, and S = 14°) [274]. The smallest probabilities 
for conversion are found, not surprisingly, around the direction parallel to the local geomagnetic field. 
Note that this conversion-free region shrinks rapidly with increasing primary energy. A similar evalua- 
tion for the Southern Site of the Pierre Auger Observatory (\B\ = 0.25 G, 1 = —35°, and S = 86°) can 
be found in [275]. 

The muonic component of an EAS differs from the electromagnetic component for two main 
reasons. First, muons are generated through the decay of cooled (E^± < 1 TeV) charged pions, and thus 
the muon content is sensitive to the initial baryonic content of the primary particle. Furthermore, since 
there is no "muonic cascade", the number of muons reaching the ground is much smaller than the number 
of electrons. Specifically, there are about 5 x 10 8 muons above 10 MeV at ground level for a vertical 
10 11 GeV proton induced shower. Second, the muon has a much smaller cross section for radiation and 




64 



E=10 195 [eV] E=10 200 [eV] E=10 205 [eV] E=10 210 [eV] 




50 100% 50 100% o 50 100% 50 100% 

Conversion Probability Conversion Probability Conversion Probability Conversion Probability 



Fig. 29: Maps of gamma ray conversion probability in the geomagnetic field for several primary energies. Az- 
imuths are as labeled, "N" denotes true north. The inner circles correspond to zenith angles 30°, 60° and horizon. 
Dashed curves indicate the opening angles of 30°, 60° and 90° to the local magnetic field. From Ref. [274]. 



pair production than the electron, and so the muonic component of an EAS develops differently than 
does the electromagnetic component. The smaller multiple scattering suffered by muons leads to earlier 
arrival times at the ground for muons than for the electromagnetic component. 

The ratio of electrons to muons depends strongly on the distance from the core; for example, 
the e + e~ to ratio for a 10 11 GeV vertical proton shower varies from 17 to 1 at 200 m from the 

core to 1 to 1 at 2000 m. The ratio between the electromagnetic and muonic shower components behaves 
somewhat differently in the case of inclined showers. For zenith angles greater than 60°, the e + e~ I /x + \±~ 
ratio remains roughly constant at a given distance from the core. As the zenith angle grows beyond 
60°, this ratio decreases, until at 9 = 75°, it is 400 times smaller than for a vertical shower. Another 
difference between inclined and vertical showers is that the average muon energy at ground level changes 
dramatically. For horizontal showers, the lower energy muons are filtered out by a combination of energy 
loss mechanisms and the finite muon lifetime: for vertical showers, the average muon energy is 1 GeV, 
while for horizontal showers it is about 2 orders of magnitude greater. The muon densities obtained in 
shower simulations using SIBYLL 2.1 fall more rapidly with lateral distance to the shower core than those 
obtained using QGSJET 01. This can be understood as a manifestation of the enhanced leading particle 
effect in SIBYLL, which can be traced to the relative hardness of the electromagnetic form factor profile 
function. The curvature of the distribution ((Pp^/dr 2 ) is measurably different in the two cases, and, 
with sufficient statistics, could possibly serve as a discriminator between hadronic interaction models, 
provided the primary species can be determined from some independent observable(s) [276]. 

High energy muons lose energy through pair production, muon-nucleus interaction, bremsstrahlung, 
and knock-on electron (5-ray) production [277]. The first three processes are discrete in the sense that 
they are characterized by high inelasticity and a large mean free path. On the other hand, because of its 
short mean free path and its small inelasticity, knock-on electron production can be considered a con- 
tinuous process. The muon bremsstrahlung cross section is suppressed by a factor of {m e /m^) 2 with 
respect to electron bremsstrahlung; see Eq. (178). Since the radiation length for air is about 36.7 g/cm 2 , 
and the vertical atmospheric depth is 1, 000 g/cm 2 , muon bremsstrahlung is of negligible importance 
for vertical air shower development. Energy loss due to muon-nucleus interactions is somewhat smaller 
than muon bremsstrahlung. As can be seen in Fig. 30, energy loss by pair production is slightly more 
important than bremsstrahlung at about 1 GeV, and becomes increasingly dominant with energy. Finally, 
knock-on electrons have a very small mean free path (see Fig. 30), but also a very small inelasticity, so 
that this contribution to the energy loss is comparable to that from the hard processes. 

3.3 Paper-and-pencil air shower modeling 

Most of the general features of an electromagnetic cascade can be understood in terms of the toy model 
due to Heitler [278]. In this model, the shower is imagined to develop exclusively via bremsstrahlung 
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and pair production, each of which results in the conversion of one particle into two (see Fig. 31). As 
was previously discussed, these physical processes are characterized by an interaction length X EM w 
37.6 g/cm 2 . One can thus imagine the shower as a particle tree with branches that bifurcate every X EM , 
until they fall below a critical energy, eo ~ 86 MeV, at which point energy loss processes dominate. 
Up to eo, the number of particles grows geometrically, so that after n = X/X EM branchings, the total 
number of particles in the shower is N w 2 n . At the depth of shower maximum X max , all particles are 
at the critical energy, eo, and the energy of the primary particle, Eq, is split among all the Af max = Eq/co 
particles. Putting this together, we get: 



ln(£ /e ) 
In 2 



^max ~ ^em — ttt; — ■ (184) 



Changes in the mean mass composition of the CR flux as a function of energy will manifest as changes 
in the mean values of X max . This change of X max with energy is commonly known as the elongation 
rate [279]: 

For purely electromagnetic showers, X max (£?) « X EM \n(E/eo), and hence D e « Xem- For con- 
venience, the elongation rate is often written in terms of energy decades, Dio = <9(X max ) / 'dlogE, 
where Dio = 2.3D e . The elongation rate obtain from the Heitler model, Dj « 84 g/cm 2 , is in very 
good agreement with the results from Monte Carlo simulations. However, the prediction for the particle 
number at maximum is overestimated by a factor of about 2 to 3. Moreover, Heitler' s model predicts a 
ratio of electron to photons of 2, whereas simulations and direct cascade measurements in the air show 
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Fig. 31: Schematic views of (a) an electromagnetic cascade and (b) a hadronic shower. In the hadron shower, 
dashed lines indicate neutral pions which do not re-interact, but quickly decay, yielding electromagnetic subshow- 
ers (not shown). Not all pions lines are shown after the n = 2 level. Neither diagram is to scale. From Ref. [280]. 



a ratio of the order of 1/6. This difference is due to the fact that multiple photons are emitted during 
bremsstrahlung and that electrons lose energy much faster than photons do. 

As we have seen, baryon-induced showers are also dominated by electromagnetic processes, thus 
Heitler's toy model is still enlightening for such cases. For proton primaries, the multiplicity rises with 
energy, and the resulting elongation rate becomes smaller. This can be understood by noting that, on 
average, the first interaction is determined by the proton mean free path in the atmosphere, A p _ a i r = Xq. 
In this first interaction the incoming proton splits into (n(E)) secondary particles, each carrying an 
average energy E/{n(E)). Assuming that X max (E) depends dominantly on the first generation of 7 
subshowers, the depth of maximum is obtained as in Eq. (184), 



X max (£) « X + X EM \u[E/{n(E))] 



(186) 



For a proper evaluation of X max , it would be necessary to sum each generation of subshowers carefully 
from their respective points of origin, accounting for their attenuation near and after the maxima. If we 
now further assume a multiplicity dependence (n(E)) « noE A , then the elongation rate becomes, 
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which corresponds to the form given by Linsley and Watson [31], 
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(188) 



(189) 



A precise calculation of a proton shower evolution has been carried out by Matthews [280], using 
the simplifying assumption that hadronic interactions produce exclusively pions, 2N W charged and 
neutral (see Fig. 31). 7r°'s decay immediately and feed the electromagnetic component of the shower, 
whereas ir^ 's soldier on. The hadronic shower continues to grow, feeding the electromagnetic component 
at each interaction, until charged pions reach a characterictic energy at which decay is more likely than 
a new interaction. The interaction length and the pion multiplicity (3N W ) are energy independent in the 
Heitler-Matthews model. The energy is equally shared by the secondary pions. For pion energy between 
1 GeV and 10 TeV, a charged multiplicity of 10 (N w = 5) is an appropriate number. 
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The first interaction diverts 1/3 of the available energy (Eq/3) into the EM component while the 
remaining 2/3 continue as hadrons. The number of hadrons increases through subsequent generation of 
particles and in each generation about 30% of the energy is transferred to the EM cascade. Therefore 
the longer it takes for pions to reach the characteristic energy £J ~ 20 GeV (below which they will 
decay into muons), the larger will be the EM component. Consequently, in long developing showers the 
energy of the muons from decaying pions will be smaller. In addition, because of the density profile of 
the atmosphere, ^ is larger high above ground than at sea level and deep showers will produce fewer 
muons. 

This positive correlation introduces a link between the primary cosmic ray interaction cross section 
on air and the muon content at ground level. According to those principles, primaries with higher cross 
sections will have a larger muon to electron ratio at ground level. 

To obtain the number of muons in the shower, one simply assumes that all pions decay into muons 
when they reach the critical energy: = (2N 7T ) Uc , where n c = \u(Eq/^) / ln(3A^-) is the number of 
steps needed for the pions to reach £J . Introducing (3 = ln(2A^-) / ln(3A^-) we have 

N» = (E /g)P . (190) 

For N n = 5, (3 = 0.85. Unlike the electron number, the muon multiplicity does not grow linearly with 
the primary energy, but at a slower rate. The precise value of (3 depends on the average pion multiplicity 
used. It also depends on the inelasticity of the hadronic interactions. Assuming that only half of the 
available energy goes into the pions at each step (rather than all of it, as done above) would lead to 
(3 = 0.93. Detailed simulations give values of f3 in the range 0.9 to 0.95 [145]. 

The first interaction yields N 1 = 2A^.o = N n ± . Each photon initiates an EM shower of energy 
Eo/(3N 7r ±) = Eq/(6N 7T ). Using pp data [15], we parametrized the charged particle production in the 
first interaction as N n ± = 41.2(£'o/l PeV) 1 / 5 . Now, from the approximation in (186), based on the sole 
evolution of the EM cascade initiated by the first interaction, we obtain 

XL* = Xo + X^HEo/iGN^eo)} 

= (470 + 581ogio[£;o/lPeV])g/cm 2 . (191) 

This falls short of the full simulation value by about 100 g/cm 2 [280]. 



EXERCISE 3.2 The depth of shower maximum obtained in Eq. (191) is only approximate since 
it considers just the first interaction as hadronic in nature. Extend the approximation to include hadronic 
interactions in the second generation of particles. Try the generalization to include all generations of 
hadronic collisions until the charged pions cool down below the critical energy. 



A good approximation of the elongation rate can be obtained when introducing the cross section 
and multiplicity energy dependence. Using a p-air cross section of 550 mb at 10 9 GeV and a rate of 
change of about 50 mb per decade of energy leads to [281] 

X ~ 90 - 9 log (Eo/EeV) g/cm 2 . (192) 

Now, assuming (as in [280]) that the first interaction initiates 2N n EM cascades of energy Eq/QN^, with 
N n oc (^o/PeV) 1 / 5 for the evolution of the first interaction multiplicity with energy, we can calculate 
the elongation rate 



yp ^ ^max = d(X In2 + X EM ln[£ /(6A^e )] 
10 dlogE dlog^o 



jjP — -iiidA -y— u — ~ ■ -EM L u / V— (193) 

or 
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Dl Q = ^Dj - 9 In 2 - 62 g/cm 2 . (194) 
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This result is quite robust as it only depends on the cross section and multiplicity evolution with energy. 
It is in good agreement with Monte Carlo simulation [145]. 

To extend this discussion to heavy nuclei, we can apply the superposition principle as a reasonable 
first approximation. In this approximation, we pretend that the nucleus comprises unbound nucleons, 
such that the point of first interaction of one nucleon is independent of all the others. Specifically, a 
shower produced by a nucleus with energy E A and mass A is modeled by a collection of A proton 
showers, each with A -1 of the nucleus energy. Modifying Eq. (184) accordingly one easily obtains 
I max oc \u(Eq/A). Assuming that B is not changing with energy, one obtains for mixed primary 
composition [31] 

d(]nA)' 



D e = X (1 - B) 



1 - 



d\nE 



(195) 



Thus, the elongation rate provides a measurement of the change of the mean logarithmic mass with 
energy. One caveat of the procedure discussed above is that Eq. (188) accounts for the energy dependence 
of the cross section and violation of Feynman scaling only for the first interaction. Note that subsequent 
interactions are assumed to be characterized by Feynman scaling and constant interaction cross sections; 
see Eq. (189). Above 10 7 GeV, these secondary interactions play a more important role, and thus the 
predictions of Eq. (195), depending on the hadronic interaction model assumed, may vary by up to 
20% [145]. In Fig. 6 we show the variation of (X max ) with primary energy as measured by HiRes and 
Auger together with predictions from a variety of hadronic interaction models. 

The muon content of an EAS at ground level 7V M , as well as the ratio N^/N e , are sensitive to 
primary composition (here, N e is the electron content at ground level). To estimate the ratio of the 
muon content of nucleus-induced to proton-induced showers, we can resort again to the principle of 
superposition. Using f3 = 0.93 we find that the total number of muons produced by the superposition of 
A individual proton showers is, oc A(E A /A) - 93 . Consequently, in a vertical shower, one expects 
a cosmic ray nucleus to produce about A om more muons than a proton. This implies that a shower 
initiated by an iron nucleus produces about 30% more muons than a proton shower. Note, however, 
that a change in the hadronic interaction model could produce a much larger effect than a change in the 
primary species. For example, replacing QGSJET 01 with SIBYLL 1.6 as the hadronic interaction model 
leads to a prediction of 60% more muons for an iron shower than for a proton shower [282]. 

The situation for gamma-induced showers is a bit different. In this case the muon component of 
the shower does not simply follow Eq. (190) because of the LPM and geomagnetic field effects [283]. 
Competition between the two processes leads to a complex behavior in N^/Nfl, as shown in Fig. 32. 

While these toys models are very useful for imparting a first intuition regarding global shower 
properties, the details of shower evolution are far too complex to be fully described by a simple analytical 
model. Full Monte Carlo simulation of interaction and transport of each individual particle is required 
for precise modeling of the shower development. At present, two Monte Carlo packages are available 
to simulate EASs: CORSIKA (COsmic Ray Simulation for KAscade) [284] and AIRES (AIR shower 
Extended Simulation) [273]. Both programs provide fully 4-dimensional simulations of the air showers 
initiated by protons, photons, and nuclei. To simulate hadronic physics, the programs make use of the 
event generators described in Sec. 3.1. Propagation of particles takes into account the Earth's curvature 
and geomagnetic field. For further details on these codes, the reader is referred to [285]. 



4 Searches for new physics beyond the electro weak scale at a/s ^ 250 TeV 
4.1 General idea 

If new physics interations occur at LHC energies, then CR collisions with cm. energies ranging up to 
250 TeV would obviously involve new physics as well. The question is, can new physics be detected by 
CR experiments? At ultrahigh energies, the cosmic ray luminosity ~7x 10 -10 (E/PeV)~ 2 cm -2 s _1 
(taking a single nucleon in the atmosphere as a target and integrating over 2tt sr) is about 50 orders of 
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Fig. 32: Ratio of the muon content for EASs produced by primary gammas and protons. The geomagnetic field is 
set to the PAO Southern site. From Ref. [283]. 



magnitude smaller than the LHC luminosity. This renders the hunt for physics beyond the electroweak 
scale futile in hadronic cosmic ray interactions occuring in the atmosphere. However, there is still a 
possibility of uncovering new physics at sub-fermi distances in cosmic neutrino interactions. 

Neutrinos are unique probes of new physics, as their interactions are uncluttered by the strong and 
electromagnetic forces and, upon arrival at the Earth, they may experience collisions with cm. energies 
up to < 250 TeV. However, rates for new physics processes are difficult to test since the flux of cosmic 
neutrinos is virtually unknown. Interestingly, it is possible in principle to disentangle the unknown flux 
and new physics processes by using multiple observables [286,287]. 

For example, possible deviations of the neutrino-nucleon cross section due to new non-perturbative 
interactions 9 can be uncovered at Auger by combining information from Earth- skimming and quasi- 
horizontal showers. In particular, if an anomalously large rate is found for deeply developing quasi- 
horizontal showers, it may be ascribed either to an enhancement of the incoming neutrino flux, or 
an enhancement in the neutrino-nucleon cross section (assuming non-neutrino final states dominate). 
However, these possibilities can be distinguished by comparing the rates of Earth- skimming and quasi- 
horizontal events. For instance, an enhanced flux will increase both quasi-horizontal and Earth- skimming 

9 Throughout these lectures we use this term to describe neutrino interactions in which the final state energy is dominated 
by the hadronic component. We are not considering here new "perturbative" physics e.g. (softly broken) supersymmetry at the 
TeV scale which would have quite different signatures in cosmic ray showers. 
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event rates, whereas an enhanced interaction cross-section will also increase the former but suppress the 
latter, because the hadronic decay products cannot escape the Earth's crust. Essentially this approach 
constitutes a straightforward counting experiment, as the detailed shower properties are not employed to 
search for the hypothesized new physics. 

The question we would like to answer is then how many Earth- skimming and quasi-horizontal 
events would we need to observe at Auger to make a convincing case for the existence of non-perturbative 
physics in which the final state is dominated by hadrons. The analysis techniques described herein consti- 
tute an entirely general approach to searching for non-perturbative interactions without any dependence 
on what hypothetical mechanism might actually cause the "hadrophilia." In Appendix D we illustrate 
one possible new physics process which may be accessible using these techniques at Auger, and which 
is beyond the reach of the LHC. 

4.2 v acceptance and systematic uncertainties 

Detailed Monte Carlo simulations are used to compute the acceptance for Earth- skimming (ES) and 
quasi-horizontal (QH) events. Neutrinos are propagated through the atmosphere, the Earth's crust, and 
the Andes mountains using an extended version [288] of the code ANIS [289]. In the simulations, the 
vN cross-sections from reference [290] are employed. Particles resulting from vN interactions are fed 
to PYTHIA [291] and r decays are simulated using TAUOLA [292]. 

The flux, energy and decay vertices of outgoing leptons are calculated inside an "active detector" 
volume of 3, 000 x 10 km 3 , including the real shape of the surface array. A relief map of the Andes 
mountains was constructed using digital elevation data from the Consortium for Spatial Information 
(CGIAR-CSI) [293]. The map of the area around the Auger site is depicted in Fig. 33. 




Fig. 33: Topography in the vicinity of the Auger site. The surface array is centered at X = Y = 0. From Ref. [294]. 

To study the response of the detector, the outputs of PYTHIA and/or TAUOLA are used as input 
for the AIRES [295] air shower simulation package. The response of the surface detector array is sim- 
ulated in detail using the Auger Offline simulation package [296]. Atmospheric background muons are 
also simulated in order to study the impact on neutrino identification, as such accidental muons can be 
wrongly classified as shower particles. The background from hadronic showers above 10 8 GeV is esti- 
mated to be 0(1) in 20 years, so for the energy bin considered in this analysis, 9.5 < log 10 (^/GeV), 
the background is negligible. 



71 



To establish benchmark neutrino rates, we use the Waxman-Bahcall bound [144] for the flux, 
<3>wb = 2.3 x 10 -8 E~ 2 GeV -1 s _1 cm -2 sr _1 , and employ the acceptance computed by the simulations 
described above. In order to estimate the systematic uncertainty associated with our lack of knowledge 
of the dependence of the flux on energy, we consider several scenarios which plausibly bracket the range 
of possibilities: 

(C/E )E-\ 
C Eq e~ 3 , 

CE~ 2 exp[- log 10 (E u /E ) 2 /(2a 2 )}, 

where C = 2.3 x 1CT 8 GeV -1 s -1 cm- 2 sr" 1 , E = 10 10 GeV, a = 0.5 GeV. The expected rates for the 
entire range over which Auger is sensitive are given in Table 6 and the rates for the high energy bin are 
given in Table 7. 

Table 6: Expected events per year at Auger in the energy range 8 < log 10 (£^/GeV), for various incident 
zenith angle (0) ranges, assuming the Waxman-Bahcall flux. 



flux 


up-going 


down-going 


ratio 


(2) 
(2) 


9 N Ut 


9 N Ue N Ut N V/l N VM 




90-95 0.68 
90-95 0.68 


60-90 0.134 0.109 0.019 0.262 
75-90 0.075 0.071 0.011 0.157 


2.58 
4.27 



Table 7: Expected events per year (Ni) at Auger in the energy range 9.5 < log 10 (^/GeV) < 10.5, for various 
incident zenith angle (0) ranges and the 4 flux models considered. 



flux 


up-going 


down-going 


ratio 




9 




9 












(1) 


90-95 


0.14 


60-90 


0.059 


0.049 


0.011 


0.12 


1.14 


(2) 


90-95 


0.15 


60-90 


0.059 


0.049 


0.096 


0.11 


1.33 


(3) 


90-95 


0.23 


60-90 


0.079 


0.062 


0.0123 


0.15 


1.53 


(4) 


90-95 


0.12 


60-90 


0.046 


0.037 


0.0080 


0.091 


1.33 


(1) 


90-95 


0.14 


75-90 


0.027 


0.031 


0.0056 


0.064 


2.14 


(2) 


90-95 


0.15 


75-90 


0.026 


0.029 


0.0048 


0.060 


2.47 


(3) 


90-95 


0.23 


75-90 


0.036 


0.041 


0.0062 


0.083 


2.75 


(4) 


90-95 


0.12 


75-90 


0.021 


0.024 


0.0040 


0.049 


2.45 



Table 8 contains a summary of systematic uncertainties on the ratio of the number of ES to QH 
events. The uncertainty in spectrum shape is taken from Table 7. The uncertainty on the parton structure 
of the nucleon is estimated by considering different PDFs (GRV92NLO [297] and CTEQ66c [298]). 
Finally, the uncertainty on the energy loss, /3 r , of r leptons as they propagate through the Earth's crust is 
derived from [299-302]. 

4.3 Auger discovery reach 

Consider first a flux of Earth- skimming r neutrinos with energy in the range 10 9 5 GeV < E v < 
10 10 - 5 GeV. The neutrinos can convert to r leptons in the Earth via the charged current interaction 



1. ^ B (E U ) = 

2. $% B (E U ) = 

3. ^ B (E U ) = 

4. ^ B (E U ) = 
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Table 8: Contributions to the systematic uncertainty on the Earth- skimming to quasi-horizontal event ratio. We 
have considered the energy range 9.5 < log 10 (E u / GeV) < 10.5 and the zenith angle range 75° < 6 < 90°. 



ratio 


flux 


PDF 




sum 




+11% 


0% 


+24% 


+ 26% 


2.47 








2.47 




-13% 


-21% 


-25% 


-35% 



v T ± N — >► r^X. In the (perturbative) SM, the interaction path length for the neutrino is 

L£ c = [NAPsa^c}- 1 , (196) 

where a^ c is the charged current cross section for a neutrino energy E v = Eq. The density of the 
material through which the neutrinos pass, p s , is about 2.65 g/cm 3 for the Earth's crust. Here we 
have neglected neutral current interactions, which at these energies only reduce the neutrino energy 
by approximately 20%, which is within the systematic uncertainty. For Eo ~ 10 10 GeV, L£ c ~ 
(9(100) km. Let us assume some hypothetical non-perturbative physics process enhances the vN cross 
section. Then the interaction path length becomes 

L t l ' ot = [iV A p s (^ c + ^p)]- 1 , (197) 

where a^ p is the non-perturbative contribution to the cross section for E v = Eq. 

Once a r is produced by a CC interaction, it can be absorbed in the Earth or escape and possibly 
decay, generating a detectable air shower. At these high energies, the r propagation length in the Earth 
is dominated by energy loss rather than the finite r lifetime. The energy loss can be expressed as 

^ = -(a T + p T E T )p 8 , (198) 

where a T characterizes energy loss due to ionization and f3 T characterizes losses through bremsstrahlung, 
pair production and hadronic interactions. At these energies, energy losses due to ionization turn out to 
be negligible, while f3 T ~ 0.8 x 10 -6 cm 2 /g [303]. From Eq. (198), we observe that the maximum path 
length for a detectable r can be written 

L T = ^ln(£ max /£ min ) , (199) 

PtPs 

where £ max « Eo is the energy at which the r is created, and i? m in is the minimal energy at which a r 
can produce a shower big enough to be detected. For ^max/^min = 10, L T = 11 km. 

The probability for a neutrino with incident nadir angle 9 to emerge as a detectable r is 

pi j 7 

P W = / T^ e ~ Z ' LXot 9 [* " ( £ ~ LT )1 ' ( 20 °) 
Jo L cc 

where £ = 2i? e cos 9 is the chord length of the intersection of the neutrino's trajectory with the Earth, 
with i?0 « 6371 km the Earth's radius. Note that we have neglected the possibility that non-perturbative 
processes could lead to a detectable signal, since the hadrons which dominate the final state will be 
absorbed in the Earth. The step function in Eq. (200) reflects the fact that a r will only escape the Earth 
if z + L T > £, as illustrated in Fig 34. 
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Fig. 34: The chord length of the intersection of a neutrino with the Earth is i = 2R§ cos 0. In the figure, the 
neutrino produces a lepton / after traveling some distance z inside the Earth's crust. If z + L r > £, the lepton will 
escape the Earth and can generate an air shower. From Ref. [179]. 



Assuming an isotropic tau neutrino flux, the number of taus that emerge from the Earth with 
sufficient energy to be detected is proportional to an "effective solid angle" 



fi e ff = J P{9) cos 9 d cos 9 d(\. 
Evaluation of the integrals [286] yields the unfortunate expression 



fLff — 27T- 



T v 



3^ T /^ot _ 1 



2R(n 




(201) 



(202) 



At the relevant energies, however, the neutrino interaction length satisfies L\ ot <C R®. In addition, if the 
hypothesized non-perturbative cross section enhancement is less than typical hadronic cross sections, we 
have L" ot » L T . With these approximations, Eq. (202) simplifies to [287] 



^eff 



tv2tt 

2n LtotL 



4i?|L£ c 



(203) 



Equation (203) describes the functional dependence of the Earth- skimming event rate on the non- 
perturbative cross section. This rate is, of course, also proportional to the neutrino flux $^ a11 at Eq. Thus, 
the number of Earth- skimming neutrinos is given by 



Nes ~ C E s 



Cj^all 



T CC 



(204) 



{(T cc -h cr Np j 

where Ces is the number of Earth- skimming events expected for some benchmark flux <3>Q a11 in the 
absence of new physics. 

In contrast to Eq. (204), the rate for quasi-horizontal showers has the form 



Nqh = Cqh 



(205) 
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Fig. 35: Event rates for Earth- skimming (left) and quasi-horizontal (right) events in the $^ all /$Q a11 — ctnp/o"cc 
plane. Note that the contours are roughly orthogonal, and so the two types of event provide complementary 
information about flux and cross section. 




Fig. 36: Projected determination of neutrino fluxes and cross sections at y/s ^ 250 TeV from future Auger data. 
The different shaded regions indicate the 90%, 95%, 99% and 3a confidence level contours in the $^ all /$Q a11 — 
ctnp/ctcc plane, for 7V°jf = 1, N$g = 10 (left), N°f = 1, N$g = 7 (middle), and N°f = 1, = 5 (right). 
The dashed line indicates the result of including the systematic uncertainty on the NLO QCD CC neutrino-nucleon 
cross section. FromRef. [294]. 



where Cqh is the number of quasi-horizontal events expected for flux $Q a11 . 

Now, without loss of generality, we normalize the neutrino flux to the Waxmann-Bahcall bound, 
i.e. <3>Q a11 = ^wb- We use ^ e ex P ecte d rates f° r the the benchmark flux to determine the values of 
Ces and Cqh in Eqs. (204) and (205), (Ces = 0.15 and Cqh = 0.06, as shown in Table 7). Given a 
flux $^ a11 and new non-perturbative physics cross section a^ p , both TVes and A/qh are determined. On 
the other hand, given just a quasi-horizontal event rate A/qh, it is impossible to differentiate between an 
enhancement of the cross section due to non-perturbative physics and an increase of the flux. However, in 
the region where significant event rates are expected, the contours of A/qh and N^s, given by Eqs. (204) 
and (205), are more or less orthogonal and provide complementary information. This is illustrated in 
Fig. 35. With measurements of Nq^ and N^ s , both and $^ a11 may be determined independently, 
and neutrino interactions beyond the (perturbative) SM may be unambiguously identified. 



75 



We now turn to determining the projected sensitivity of Auger to neutrino fluxes and cross sections. 
The quantities TVes and TVqh as defined in Eqs. (204) and (205) can be regarded as the theoretical values 
of these events, corresponding to different points in the cfr^aii/cfr^aii _ a^p/vcc parameter space. For a 
given set of observed rates N^ s and Nq^, two curves are obtained in the two-dimensional parameter 
space by setting A^jf = A^es and Nq^ = Aqh- These curves intersect at a point, yielding the most 
probable values of flux and cross section for the given observations. Fluctuations about this point define 
contours of constant x 2 in an approximation to a multi-Poisson likelihood analysis [304]. The contours 
are defined by 



where i = ES, QH [305]. In Fig. 36, we show results for three representative cases. Assuming (iVgg s = 1, 
Aqjj = 10), (N°f = 1, iVgg = 7)? and (7V obs = 1? jyobs = 5) we show the 90%9 95%9 99% and 3a C L 

contours for 2 d.o.f. (x 2 = 4.61, 5.99, 9.21, and 11.83, respectively). For N°f = 1 and A^ s = 10, the 
possibility of a SM interpretation along the a^ p = axis (taking into account systematic uncertainties) 
would be excluded at greater than 99% CL for any assumed flux. 

In summary, we have found that observation of 1 Earth- skimming and 10 quasi-horizontal events 
would exclude the (perturbative) SM at the 99% CL. Thus the expected low neutrino "luminosity" is not 
at all a show-stopper, and the observatory has the potential to uncover new physics at scales exceeding 
those accessible to the LHC. If new non-perturbative physics exists, a decade or so would be required to 
uncover it at Auger in the best case scenario (cosmic neutrino flux at the Waxman-Bahcall level and vN 
cross section about an order of magnitude above the SM prediction). Of course future CR experiments 
should benefit from a much larger aperture, making such discovery conceivable in a short time scale. An 
optimist may even imagine the possibility of probing QCD dynamics [306]. 
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Appendices 

A Cosmogenic /3-DK and A* processes 

If UHECRs are heavy nuclei, the relic photons can excite the GDR at nucleus energies E > 10 11 GeV 
(corresponding to 10 MeV — 30 MeV in the nuclear rest frame), and thus there should be 

accompanying photo-dissociated free neutrons, themselves a source of /3-decay antineutrinos. The decay 
mean free path of a neutron is cY n r n = (E n /10 u GeV) Mpc, the lifetime being boosted from its rest 
frame value r n = 886 s to its lab value via T n = E n /m n . Compared to cosmic distances > 100 Mpc, 
the decay of even the boosted neutron may be taken as nearly instantaneous, and thus all free neutrons are 
themselves a source of /3-decay cosmogenic antineutrinos. The neutron emissivity Q n (E n ), defined as 
the mean number of neutrons emitted per comoving volume per unit time per unit energy as measured at 
the source can be estimated as follows. Neutrons with energies above 10 9 3 GeV have parent iron nuclei 
with r > Tq ~ 2 x 10 9 which are almost completely disintegrated in distances of less than 100 Mpc (see 




(206) 
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Sec. 2.3.2). Thus, it is reasonable to define a characteristic time r r given by the moment at which the 
number of nucleons is reduced to 1 je of its initial value A, and presume the nucleus, emitted at distance d 
from the Earth, is a traveling source that at D « (d — cr r ) disintegrates into A nucleons all at once [120]. 
Then, the number of neutrons with energy E n = Ea/A can be approximated by the product of 1/2 the 
number of nucleons generated per nucleus and the number of nuclei emitted, i.e. Q n (E n ) =NQa, where 
N = A — Z is the mean neutron number of the source nucleus. Now, to obtain an estimate of the 
diffuse antineutrino flux one needs to integrate over the population of nucleus-emitting-sources out to 
the horizon [307, 308] 



1 f 1T , ^ ^ ^ \ ( Dm n \] f Q 7 dP , N 

tio J |_ \ E n r n J ] J dev 

6 [E^- E n e 17 (l + cos O^/rrin] , (A.l) 



47T 

^ dcosOjj 



where the r 2 in the volume integral is compensated by the usual 1/r 2 fall-off of flux per source. Here, 
Hq is the Hubble constant, E v and E n are the antineutrino and neutron energies in the lab, 6 V is the 
antineutrino angle with respect to the direction of the neutron momentum in the neutron rest-frame, 
and e„ is the antineutrino energy in the neutron rest-frame. The last three variables are not observed 
by a laboratory neutrino-detector, and so are integrated over. The observable E v is held fixed. The 
delta- function relates the neutrino energy in the lab to the three integration variables. The parameters 
appearing in Eq. (A.l) are the neutron mass and rest-frame lifetime (m n and r n ). Finally, dP/dejj 
is the normalized probability that the decaying neutron produces a V with energy in the neutron 
rest-frame. Note that the maximum V energy in the neutron rest frame is very nearly Q = m n — 
m p — m e = 0.71 MeV. Integration of Eq. (A.l) can be easily accomplished, especially when two good 
approximations are applied [309] . The first approximation is to think of the /3-decay as a 1 —> 2 
process of 5m n -> e~ + 17, in which the neutrino is produced monoenergetically in the rest frame, with 
ej7 = eo ~ SrriN(l — m 2 e / 5 2 rriN) /2 ~ 0.55 MeV, where Stun ~ 1.30 MeV is the neutron-proton mass 
difference. In the lab, the ratio of the maximum V energy to the neutron energy is 2eo/m n ~ 10 -3 , 
and so the boosted V has a spectrum with E v e (0, 10 -3 E n ). The second approximation is to replace 
the neutron decay probability 1 - e - Dm n/E n r n w ^ a ste p f unc ti n 6(£^f ax - E n ) at some energy 
^max ^ q^jj mn /r n ) = (D/10 Mpc) x 10 12 GeV Combining these two approximations we obtain 

where E™ m = max{^ ; r , A ^ Ev }, and £^™ ax is the energy cutoff at the nucleus-emitting-source 
< A(D/W Mpc) x 10 12 GeV. For Q A oc E^ a , integration of Eq. (A.2) leads to 

Ea,f 



^ e (Eu) w 10 6 
where Ep > 10 6 3 (56/A) GeV and 



ZTTnax 
A 



rmax \ csl 



The sub-PeV antineutrino spectrum is flat, as all the free neutrons have sufficient energy E n > Er /A, 
to contribute equally to all the V energy bins below 10 6 GeV. Taking a = 2 as a reasonable example, 

and inputting the observational value for iron nuclei, E\ ro Jqr 

becomes 



w 10 5 eV m~ 2 s" 1 sr" 1 , Eq. (A.3) 

ro 



El <& Ve w 4 x 10 1 ( ^ ) eVm" 2 s" 1 sr 1 . (A.5) 
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Note that the /3-decay process gives initial antineutrino flavor ratios 1:0:0 and Earthly ratios nearly 3 : 
1:1. Compared to full-blown Monte Carlo simulations [310,311], this back of the envelope calculation 
underestimates the flux by about 30%. This is because the preceeding calculation does not account for 
possible neutrons produced in £>7cmb collisions. Of course the situation described above represents the 
most extreme case, in which all cosmic rays at the end of the spectrum are heavy nuclei. A more realistic 
guess would be that the composition at the end of the spectrum is mixed. 

Photodisintegration of high energy nuclei is followed by immediate photo-emission from the ex- 
cited daughter nuclei. For brevity, we label the photonuclear process A + 7 — > A'* + X, followed by 
A f * —> A' + 7-ray as "A*" [312]. The GDR decays dominantly by the statistical emission of a single 
nucleon, leaving an excited daughter nucleus (A — 1)*. The excited daughter nuclei typically de-excite 
by emitting one or more photons of energies e^ xn ^1 — 5 MeV in the nuclear rest frame. The lab-frame 
energy of the 7-ray is then ei^ AB = Ta ^ xn , where Ta = Ea/ttia is the boost factor of the nucleus in 
the lab. 

Of interest here is the 7-ray flux produced when the photo-dissociated nuclear fragments produced 
on the CMB and CIB de-excite. These 7-rays create chains of electromagnetic cascades on the CMB and 
CIB, resulting in a transfer of the initial energy into the so-called Fermi-LAT region, which is bounded 
by observation to not exceed cj cas ~ 5.8 x 10 _7 eV/cm 3 [170]. Fortunately, we can finesse the details 
of the calculation by arguing in analogy to work already done. The photodisintegration chain produces 
one /3-decay neutrino with energy of order 0.5 MeV in the nuclear rest frame, for each neutron produced. 
Multiplying this result by 2 to include photodisintegration to protons in addition to neutrons correctly 
weights the number of steps in the chain. Each step produces on average one photon with energy ~ 
3 MeV in the nuclear rest frame. In comparison, about 12 times more energy is deposited into photons. 
Including the factor of 12 relating uo 1 to ojp e , we find from (A.5) that cosmogenic photodisintegration/de- 
excitation energy, cj cas ^ 1.4 x 1CT 11 (56/A) eV/cm 3 , is more than three orders of magnitude below 
the Fermi-LAT bound. 10 This result appears to be nearly invariant with respect to varying the maximum 
energy of the Fe injection spectrum (with a larger £ max , the additional energy goes into cosmogenic pion 
production). Thus, there is no constraint on a heavy nuclei cosmic-ray flux from the A* mechanism. 

Neutron emission from a cosmological distribution of CR sources would also lead to a flux of V e . 
In analogy to our previous estimate, we sum over the neutron-emitting sources out to the edge of the 
universe; Eq. (A.2) becomes [314] 



where Q n {E n ) is the neutron volume emissivity. An upper limit can be placed on Q n by assum- 
ing an extreme situation in which all the CRs escaping the source are neutrons, i.e., ej^ > 10 1 = 
J dE n E n Q n (E n ). With the production rate of ultrahigh energy protons (91), and an assumed injection 
spectrum Q n oc E~ 2 , Eq. (A.6) gives 



which is about three orders of magnitude below the Waxman-Bahcall bound. 
B Energy density of the EM cascade 

EM interactions of photons and leptons with the extragalactic radiation backgrounds and magnetic field 
can happen on time-scales much shorter than their production rates. The relevant processes with back- 
ground photons contributing to the differential interaction rates 7 ee , 7 7e , and 7 e7 are inverse Comp- 
ton scattering (ICS), e ± + 7b gr —> e ± + 7, pair production (PP), 7 + 7b gr —> e + + e~, double pair 

10 These rough calculations are consistent with the more rigurous analysis presented in [313]. 




(A.6) 



El ^ w 3 x 10" 11 GeVcm^s" 1 sr 



-1 



(A.7) 
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production (DPP) 7 + 7b gr —> e + + e + e + + e , and triple pair production (TPP), e ± + 7b gr —> 
e ± + e+ + e~ [106,315]. 

High energetic electrons and positrons may also lose energy via synchrotron radiation in the in- 
tergalactic magnetic field B with a random orientation sin 9 with respect to the velocity vector. We will 
assume in the following that the field strength B = 10 -12 G. 11 This leads to an efficient transfer of 
energy into the EM cascade. The synchrotron power spectrum (W eV _1 ) has the form 



P(£ e , £? 7 ) = F(Ej/E c ) ; F(t) = i / <fe tf 5/3 (z) , 

Z7T m e j t 



(B.l) 



where K 5 / 3 is the modified Bessel function and E c = (3e J Bsin0/2m e )(^ e /m e ) 2 . This can be treated 
as a continuous energy loss of the electrons and positrons with a parameter 12 



W*e) = ^ / decs* / ^ 7 W,£ 7 ) = f (^) (^) • (B.2) 

We will assume in the calculation that the intergalactic magnetic field is primordial with a (flux-conserving) 
redshift dependence B(z) = (1 + z) 2 B(0). Note that the synchrotron energy loss has then a redshift 
dependence similar to BH pair production in the CMB, i.e. b syn (z, E) = (1 + ^) 2 6 syn (0, (1 + z)E). It is 
also convenient to define 7^ n (i? e , E 7 ) = V(E ei E 1 )/E 1 , which has an analogous redshift dependence, 
i.e. ffi(z, E e , = (1 + z) 4 7 !7 n (0> (1 + z)E e , (1 + z)E 1 ). 

The fast evolution of the cascade is governed by the set of differential equations, 

d i Y e {E)=-T e (E)Y e {E) + d E {b{E)Y e {E)) + J dE' [ lie {E' , E)Y 1 {E') + lee (E' , E)Y e {E')] , 

d t Yy(E) = - T 7 (E)Y 7 (E) + J dE , V ^ E ^ Y e (E') + J dE'^ ei (E r , E)Y e (E r ) , (B.3) 

which determines the evolution on short time-scales At Y pi <C 1 (the redshift z is kept fixed meanwhile). 
The initial condition Yy/ e (E)\f =0 is given by the sum of previously developed cascades and the newly 
generated contributions from proton interactions. 

The solution of the system (B.3) for an infinitesimally small step At can be written for a discrete 
energy spectrum, Ni ~ AEiYi, as 

With the transition matrix T(At), defined by Eq. (B.4), we can efficiently follow the development of the 
EM cascade over a distance At = 2 N At via matrix doubling [167]: 

T{2 N At)~ [T(At)] 2N . (B.5) 
The total energy of the cascade can be obtained by integrating Eq. (120): 

J dE E^ s (z, E) = - J dEEdg [6 cas (z, S)2t v {z, E)\ . (B.6) 

Integrating the r.h.s. by parts yields 

r.h.s. = - J dEd E E^-b^{z,S)^ v {z,E) + J dE^-b cas (z^)^ p (z, E) . (B.7) 

n For B > 10~ 12 G, TPP by electrons can be neglected [316]. For E < 10 12 GeV, DPP of photons can also be safely 
neglect in the calculation [317]. 

12 Note the identity f dE [E fe(6n e ) + f dE'V(E\ E)n e ] = 0, implying overall energy conservation. 
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E [GeV] 

Fig. B.l: Cosmogenic photon (blue lines) and neutrino (red line) fluxes from a "vanilla" proton test- spectrum 
(shown as black line) of UHECRs, assuming two exponential cutoffs (E m i n = 10 8 GeV and E max = 10 11,5 GeV), 
power index 7 = 2 and evolution with n = 3 and < z < 2. The fluxes are computed for 3 extragalactic magnetic 
field strengths: 1 pG, 1 nG, and 50 nG. Proton diffusion in extragalactic magnetic fields B » nG has not yet been 
included in the calculation. We have neglected radiative losses of muons and pions before decaying into neutrinos; 
the energy loss in muons is down by a factor 10 6 compared to electrons. This figure is courtesy of Markus Ahlers. 



The first term vanishes since 6 cas = for sufficiently low energies and 3? v = beyond the maximal 
energy. The time integration of the l.h.s. between the present epoch (t = 0) and the first sources (t max ) 
gives 

J dt[l.h.s.] = J dEEn cas (E) = u cas , (B.8) 

hence we obtain Eq. (121). 

Synchrotron radiation in strong magnetic fields, 1 nG < B < 50 nG, can also by-pass the EM cas- 
cade and transfer energy to sub-GeV photons that are unconstrained by the Fermi-LAT spectrum [318]. 
This would be relevant for electrons around E e ~ 10 9 GeV, where synchrotron loss starts to dominate 
over ICS loss in the CMB; see Fig. B.l. From Eqs. (20) and (27) we can verify that (E* yn ) < 1 GeV. Of 
course, reaching a 50 nG field requires some astrophysical generator to augment the primordial field. 13 
For simplicity, we have taken the redshift evolution of the generator to be the same as the evolution of 
the primordial field. If were to assume an evolution which follows SFR, this will not significantly alter 
our conclusions. 

C TOTal Elastic and diffractive cross section Measurement 

The TOTEM experiment is dedicated to the measurement of the total pp cross section with the luminosity 
independent method based on the Optical Theorem, 

a tot = Im/(0) , (C.l) 



where /(#) satisfies 



at s ail s 



13 For Xb ~ 1 Mpc, the current upper limit on primordial seed fields is about 1 nG [319, 320]. 
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with i) the angle of the scattered proton with respect to the beam direction. Squaring Eq. (C.l) we obtain 



2 16tt Im 2 /(0) da el 

a i — 



'tot 



Re 2 /(0) + Im 2 /(0) dt 



16tt [dN el /dt] t=0 _ 



t=o 



1 + P 2 C 



where C is the integrated luminosity. Now, following [321, 322], we can obtain the total cross section 
independently from C, by using a tot = (a e \ + cr ine i) = (N e \ + 7V ine i)/£, 

167T [dN el /dt] t=0 

1 + p 2 (N el + 7V ine i) 

Here, N e \ and N[ ne \ are the numbers of elastic and inelastic events, and p = 0.10 ± 0.01 is the ratio 
between the real and imaginary parts of the forward scattering amplitude [323]. 14 The difficult aspect 
of this measurement is obtaining a good extrapolation of the cross section for low momentum transfer. 
Recall that —t = s(l — cos #)/2 ~ si} 2 / 4, and so t — >> implies a measurement in the extreme forward 
direction. The TOTEM experiment aims to measure down to values of \t\ ~ xlO -3 GeV 2 , which 
corresponds to i? « 4.5 /irad [324]. 



D Raiders of the lost holy grail 

In 1976 't Hooft observed that the Standard Model does not strictly conserve baryon and lepton num- 
ber [325, 326]. Rather, non-trivial fluctuations in SU{2) gauge fields generate an energy barrier interpo- 
lating between topologically distinct vacua. An index theorem describing the fermion level crossings in 
the presence of these fluctuations reveals that neither baryon nor lepton number is conserved during the 
transition, but only the combination B — L. Inclusion of the Higgs field in the calculation modifies the 
original instanton configuration [327]. An important aspect of this modification (called the "sphaleron") 
is that it provides an explicit energy scale of about 10 TeV for the height of the barrier. This barrier can 
be overcome through thermal transitions at high temperatures [328-330], providing an important input to 
any calculation of cosmological baryogenesis. More speculatively, it has been suggested [331-333] that 
the topological transition could take place in two particle collisions at very high energy. The anomalous 
electroweak contribution to the partonic process can be written as 

di{s) = 5.3 x 10 3 mb • e ~^ /aw) Fw ^ e) , (D.l) 

where aw — 1/30, the tunneling suppression exponent i^V(e) is sometimes called the "holy-grail func- 
tion", and e = >/§ / \\.i\mw I otw) — V^/30 TeV. Thus, it is even possible that at or above the sphaleron 
energy the cross-section could be of (9(mb) [334]. Of particular interest here would be an enhancement 
of the neutrino cross section over the perturbative SM estimates, say by an order of magnitude in the 
energy range 9.5 < log 10 (^/GeV) < 10.5. 

The argument for strong damping of the anomalous cross section for y/§ > 30 TeV was convinc- 
ingly demonstrated in [335, 336], in the case that the classical field providing the saddle point interpo- 
lation between initial and final scattering states is dominated by spherically symmetric configurations. 
This 0(3) symmetry allows the non- vacuum boundary conditions to be fully included in extremizing 
the effective action. In [337] it was shown that a sufficient condition for the 0(3) dominance is that the 
interpolating field takes the form of a chain of "lumps" which are well- separated, so that the each lump 
lies well into the exponentially damped region of its nearest neighbors. However, we are not aware of 
any reason that such lumped interpolating fields should dominate the effective action. It is thus of interest 
to explore the other extreme, in which non-spherically symmetric contributions dominate the effective 
action (and let experiment be the arbiter). 

14 Note that the quoted value of p is an extrapolation to = 14 TeV, and may be measured by the LHC experiments. 
Otherwise, it will contribute to the uncertainty in cr to t. 
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It was shown [334] that for the simple sphaleron configuration s-wave unitarity is violated for 
y/§ > 4ttMw/(^w ~ 36 TeV. If for y/§ > 36 TeV we saturate unitarity in each partial wave, then this 
yields a geometric parton cross section ttR 2 , where R is some average size of the classical configuration. 
As a fiducial value we take the core size of the Manton-Klinkhamer sphaleron, R ~ 4/Mw — 10 -15 cm. 
In this simplistic model, the uN cross section is 



where x min = s min /s = (36) 2 /2ME u ~ 0.065. In the region 0.065 < x < 3(0.065) the PDF for the 
up and down quarks is well approximated by / ~ 0.5 /x, so the expression for the cross section becomes 



where the last factor of 2/2 takes into account the (mostly) 2 contributing quarks (u, d) in this range of 
x, and the condition that only the left-handed ones contribute to the scattering. This is about 80 times 
the SM cross section. Of course this calculation is very approximate and the cross section can easily be 
smaller by a factor of 10 {e.g. if R is 1/3 of the fiducial value used). 

References 

[1] V. F. Hess, Phys. Z. 13, 1804 (1912). 

[2] P. Auger, R. Maze, T. Grivet-Meyer, Comptes Rendus 206, 1721 (1938). 

[3] P. Auger, P. Ehrenfest, R. Maze, J. Daudin, Robley, and A. Freon, Rev. Mod. Phys. 11, 288 (1939). 
[4] M. Nagano and A. A. Watson, Rev. Mod. Phys. 72, 689 (2000). 
[5] A. Letessier-Selvon and T. Stanev, arXiv: 1103.0031 [astro-ph.HE]. 

[6] L. Anchordoqui, T. C. Paul, S. Reucroft and J. Swain, Int. J. Mod. Phys. A 18, 2229 (2003) 

[arXiv:hep-ph/0206072]. 
[7] T. Abu-Zayyad et al. 9 Nucl. Instrum. Meth. A 450, 253 (2000). 
[8] R. M. Baltrusaitis et al., Nucl. Instrum. Meth. A 240, 410 (1985). 
[9] J. Abraham et al. [Pierre Auger Collaboration], Nucl. Instrum. Meth. A 523, 50 (2004). 
[10] J. Abraham et al. [Pierre Auger Collaboration], Nucl. Instrum. Meth. A 613, 29 (2010). 
[11] J. A. Abraham et al. [Pierre Auger Collaboration], Nucl. Instrum. Meth. A 620, 227 (2010) 
[arXiv:0907.4282]. 

[12] P. Abreu et al. [Pierre Auger Collaboration], Astropart. Phys. 34, 368 (2011) [arXiv: 1010.6162 
[astro-ph.HE]]. 

[13] Y. Takahashi et al. [JEM-EUSO Collaboration], New J. Phys. 11, 065009 (2009) [arXiv:0910.4187 

[astro-ph.HE]]. 
[14] T. Weiler, private communication. 

[15] C. Amsler et al. [Particle Data Group], Phys. Lett. B 667, 1 (2008). 
[16] M. Aglietta et al. [EAS-TOP Collaboration], Astropart. Phys. 21, 583 (2004). 
[17] T. Antoni et al. [KASCADE Collaboration], Astropart. Phys. 24, 1 (2005) [arXiv:astro- 
ph/0505413]. 

[18] W. D. Apel et al [KASCADE Collaboration], Astropart. Phys. 31, 86 (2009) [arXiv:08 12.0322 
[astro-ph]]. 

[19] J. Bluemer, R. Engel and J. R. Hoerandel, Prog. Part. Nucl. Phys. 63, 293 (2009) [arXiv:0904.0725 
[astro-ph.HE]]. 




!min partons 



(D.2) 



<lackdisk(^) - 7ri? 2 (0.5) (In 3) (2/2) 



~ 1.5 x 10" 30 cm 2 



(D.3) 



82 



[20] M. Hillas, Phys. Lett. A 24, 677 (1967). 

[21] V. Berezinsky, A. Z. Gazizov and S. I. Grigorieva, Phys. Rev. D 74, 043005 (2006) [arXiv:hep- 
ph/0204357]. 

[22] R. Aloisio, V. Berezinsky, P. Blasi, A. Gazizov, S. Grigorieva and B. Hnatyk, Astropart. Phys. 27, 

76 (2007) [arXiv:astro-ph/0608219]. 
[23] G. Cocconi, Astropart. Phys. 4, 281 (1996). 
[24] A. A. Penzias and R. W. Wilson, Astrophys. J. 142, 419 (1965). 
[25] K. Greisen, Phys. Rev. Lett. 16, 748 (1966). 

[26] G. T. Zatsepin and V. A. Kuzmin, JETP Lett. 4, 78 (1966) [Pisma Zh. Eksp. Teor. Fiz. 4, 114 
(1966)]. 

[27] J. N. Bahcall and E. Waxman, Phys. Lett. B 556, 1 (2003) [arXiv:hep-ph/0206217]. 
[28] R. U. Abbasi et al. [HiRes Collaboration], Phys. Rev. Lett. 100, 101101 (2008) [arXiv:astro- 
ph/0703099]. 

[29] J. Abraham et al. [Pierre Auger Collaboration], Phys. Rev. Lett. 101, 061101 (2008) 
[arXiv:0806.4302]. 

[30] J. Abraham et al. [Pierre Auger Collaboration], Phys. Lett. B 685, 239 (2010) [arXiv: 1002.1975 
[astro-ph.HE]]. 

[31] J. Linsley and A. A. Watson, Phys. Rev. Lett. 46, 459 (1981). 

[32] L. D. Landau and I. Pomeranchuk, Dokl. Akad. Nauk Ser. Fiz. 92, 535 (1953). 

[33] A. B. Migdal, Phys. Rev. 103, 1811 (1956). 

[34] J. Abraham et al. [Pierre Auger Collaboration], Astropart. Phys. 27, 155 (2007) [astro-ph/0606619]. 
[35] J. Abraham et al. [Pierre Auger Collaboration], Astropart. Phys. 29, 243 (2008) [arXiv:0712.1 147 
[astro-ph]]. 

[36] J. Abraham et al. [Pierre Auger Collaboration], Astropart. Phys. 31, 399 (2009) [arXiv:0903.1 127 
[astro-ph.HE]]. 

[37] R. U. Abbasi et al. [HiRes Collaboration], Phys. Rev. Lett. 104, 161101 (2010) [arXiv:09 10.4 184 
[astro-ph.HE]]. 

[38] D. R. Bergman [HiRes Collaboration], Nucl. Phys. Proc. Suppl. 136, 40 (2004) [arXiv:astro- 
ph/0407244]. 

[39] D. J. Bird et al. [HiRes Collaboration], Phys. Rev. Lett. 71, 3401 (1993). 

[40] J. Abraham et al. [Pierre Auger Collaboration], Phys. Rev. Lett. 104, 091101 (2010) 

[arXiv: 1002.0699 [astro-ph.HE]]. 
[41] R. Ulrich, R. Engel, S. Muller, F. Schussler and M. Unger, Nucl. Phys. Proc. Suppl. 196, 335 (2009) 

[arXiv:0906.3075 [astro-ph.HE]]. 
[42] L. A. Anchordoqui, M. T. Dova, L. N. Epele and S. J. Sciutto, Phys. Rev. D 59, 094003 (1999) 

[arXiv:hep-ph/9810384]. 
[43] J. Linsley, Phys. Rev. Lett. 34, 1530 (1975). 

[44] P. Sommers, Astropart. Phys. 14, 271 (2001) [arXiv:astro-ph/0004016]. 
[45] P. Sokolsky, P. Sommers and B. R. Dawson, Phys. Rept. 217, 225 (1992). 
[46] J. Aublin, E. Parizot, Astron. Astrophys. 441, 407 (2005). 
[47] R. Bonino et al. , in preparation. 

[48] P. Abreu et al. [Pierre Auger Collaboration], Astropart. Phys. 34, 267 (2011). [arXiv: 1103.2721 
[astro-ph.HE]]. 

[49] N. Hayashida et al. [AGASA Collaboration], Astropart. Phys. 10, 303 (1999) [arXiv:astro- 
ph/9807045]. 



83 



[50] T. Stanev, Astrophys. J. 479, 290 (1997) [arXiv:astro-ph/9607086]. 
[51] J. Candia, S. Mollerach and E. Roulet, JCAP 0305, 003 (2003) [arXiv:astro-ph/0302082]. 
[52] A. Calvez, A. Kusenko and S. Nagataki, Phys. Rev. Lett. 105, 091101 (2010) [arXiv: 1004.2535 
[astro-ph.HE]]. 

[53] A. H. Compton and I. A. Getting, Phys. Rev. 47, 817 (1935). 

[54] M. Kachelriess and P. D. Serpico, Phys. Lett. B 640, 225 (2006) [arXiv:astro-ph/0605462]. 
[55] J. Abraham et al. [Pierre Auger Collaboration], Science 318, 938 (2007) [arXiv:07 1 1 .2256 [astro- 
ph]]. 

[56] J. Abraham et al. [Pierre Auger Collaboration], Astropart. Phys. 29, 188 (2008) [Erratum-ibid. 30, 

45 (2008)] [arXiv:07 12.2843 [astro-ph]]. 
[57] R. U. Abbasi etal, Astropart. Phys. 30, 175 (2008) [arXiv:0804.0382 [astro-ph]]. 
[58] P. Abreu et al. [Pierre Auger Collaboration], Astropart. Phys. 34, 314 (2010) [arXiv: 1009. 1855]. 
[59] D. S. Gorbunov, P. G. Tinyakov, 1. 1. Tkachev, S. V. Troitsky, [arXiv:0804.1088 [astro-ph]]. 
[60] A. M. Hillas, Ann. Rev. Astron. Astrophys. 22, 425 (1984). 
[61] W. F. G. Swann, Phys. Rev. 43, 217 (1933). 
[62] E. Fermi, Phys. Rev. 75, 1169 (1949). 

[63] T. K. Gaisser, Cosmic Rays And Particle Physics, (Cambridge, UK: Univ. Press, 1990). 
[64] R. Blandford and D. Eichler, Phys. Rept. 154, 1 (1987). 

[65] M. Ahlers, L. A. Anchordoqui, J. K. Becker, T. K. Gaisser, F. Halzen, D. Hooper, S. R. Klein. P. 

Meszaros, S. Razzaque, and S. Sarkar, FERMILAB-FN-0847-A, YITP-SB- 10-01. 
[66] L. O. Drury, Contemp. Phys. 35, 231 (1994). 
[67] A. R. Bell, Mon. Not. Roy. Astron. Soc. 182, 147 (1978). 
[68] R. D. Blandford and J. P. Ostriker, Astrophys. J. 221, L29 (1978). 

[69] L. A. Anchordoqui, G. E. Romero and J. A. Combi, Phys. Rev. D 60, 103001 (1999) [arXiv:astro- 
ph/9903145]. 

[70] D. F. Torres and L. A. Anchordoqui, Rept. Prog. Phys. 67, 1663 (2004) [arXiv:astro-ph/0402371]. 
[71] B. L. Fannaroff and J. M. Riley, Mon. Not. Roy. Astron. Soc. 167, 31 (1974). 
[72] R. D. Blandford and M. J. Rees, Mon. Not. Roy. Astron. Soc. 169, 395 (1974). 
[73] A. Rosen, P. A. Hughes, G. C. Duncan and P. E. Hardee, Astrophys. J. 516, 729 (1999) [arXiv:astro- 
ph/9901046]. 

[74] M. C. Begelman, R. D. Blandford and M. J. Rees, Rev. Mod. Phys. 56, 255 (1984). 

[75] L. O. Drury, Rept. Prog. Phys. 46, 973 (1983). 

[76] A. N. Kolmogorov, C. R. Acad. URSS 30, 201 (1941). 

[77] J. P. Rachen and P. L. Biermann, Astron. Astrophys. 272, 161 (1993) [arXiv:astro-ph/9301010]. 
[78] G. B. Rybicki and A. P. Lightman, Radiative Processes in Astrophysics (New York: Wiley- 

Interscience, 1979). 
[79] T. A. Armstrong et al, Phys. Rev. D 5, 1640 (1972). 
[80] P. L. Biermann and P. A. Strittmatter, Astrophys. J. 322, 643 (1987). 

[81] L. A. Anchordoqui and H. Goldberg, Phys. Rev. D 65, 021302 (2002) [arXiv:hep-ph/0106217]. 
[82] F. P. Israel, Astron. Astrophys. Rev. 8, 237 (1998). 

[83] J. J. Condon, G. Helou, D. B. Sanders, and B. T. Soifer, Astrophys. J. Suppl. 103, 81 (1996). 
[84] F. Aharonian [H.E.S.S. Collaboration], Astrophys. J. 695, L40 (2009) [arXiv:0903.1582 [astro- 
ph.CO]]. 

[85] M. J. Hardcastle, D. M. Worrall, R. P. Kraft, W. R. Forman, C. Jones and S. S. Murray, Astrophys. 
J. 593, 169 (2003) [arXiv:astro-ph/0304443]. 



84 



[86] J. O. Burns, E. D. Feigelson, and E. J. Schreier, Astrophys. J. 273, 128 (1983). 

[87] G. E. Romero, J. A. Combi, L. A. Anchordoqui and S. E. Perez Bergliaffa, Astropart. Phys. 5, 279 
(1996) [arXiv:gr-qc/9511031]. 

[88] M. Honda, Astrophys. J. 706, 1517 (2009) [arXiv:09 11.0921 [astro-ph.HE]]. 

[89] N. Junkes, R. F. Haynes, J. I. Harnett, and D. L. Jauncey, Astron. Astrophys. 269, 29 (1993) [Erra- 
tum, ibid 274, 1009(1993)]. 

[90] P. Sreekumar, D. L. Bertsch, R. C. Hartman, P. L. Nolan and D. J. Thompson, Astropart. Phys. 11, 
221 (1999) [arXiv:astro-ph/9901277]. 

[91] J. E. Grindlay et al, Astrophys. J. 197, L9 (1975). 

[92] A. A. Abdo et al. [Fermi-LAT Collaboration], Science 328, 725 (2010) [arXiv: 1006.3986 [astro- 
ph.HE]]. 

[93] V. L. Ginzburg and S. I. Syrovatskii, Ann. Rev. Astron. Astrophys. 3, 297 (1965). 

[94] G. J. Fishman and C. A. Meegan, An. Rev. Astron. Astrophys. 33, 415 (1995). 

[95] B. Link and R. I. Epstein, Astrophys. J. 466, 764 (1996) [arXiv:astro-ph/9601033]. 

[96] C. A. Meegan et. al., Nature 355, 143 (1992). 

[97] M. R. Metzger et. al., Nature 387, 878 (1997). 

[98] T. Piran, Phys. Rept. 314, 575 (1999) [arXiv:astro-ph/9810256]. 

[99] P. Meszaros, Rept. Prog. Phys. 69, 2259 (2006) [arXiv:astro-ph/0605208]. 

[100] E. Waxman, Lect. Notes Phys. 598, 393 (2003). [arXiv:astro-ph/0303517]. 

[101] E. Waxman, Phys. Rev. Lett. 75, 386 (1995) [arXiv:astro-ph/9505082]. 

[102] W. Coburn and S. E. Boggs, Nature 423, 415 (2003). 

[103] E. Nakar, T. Piran and E. Waxman, JCAP 0310, 005 (2003) [arXiv:astro-ph/0307290]. 
[104] F. W. Stecker, Phys. Rev. Lett. 21, 1016 (1968). 

[105] V. S. Berezinsky and S. I. Grigor'eva, Astron. Astrophys. 199, 1 (1988). 

[106] G. R. Blumenthal, Phys. Rev. D 1, 1596 (1970). 

[107] R. M. Barnett et al. [Particle Data Group], Phys. Rev. D 54, 1 (1996). 

[108] L. Montanet et al. [Particle Data Group], Phys. Rev. D 50, 1 173 (1994). See p. 1335. 

[109] I. Golyak, Mod. Phys. Lett. A 7, 2401 (1992). 

[110] L. A. Anchordoqui, M. T. Dova, L. N. Epele and J. D. Swain, Phys. Rev. D 55, 7356 (1997) 
[arXiv:hep-ph/9704387]. 

[Ill] M. Abramowitz, I. A. Stegun, "Handbook of Mathematical Functions" (Dover, New York, 1970). 

[112] J. L. Puget, F W. Stecker and J. H. Bredekamp, Astrophys. J. 205, 638 (1976). 

[113] M. J. Chodorowski, A. A. Zdziarski, and M. Sikora, Astrophys. J. 400, 181 (1992). 

[114] S. Michalowski, D. Andrews, J. Eickmeyer, T. Gentile, N. Mistry, R. Talman and K. Ueno, Phys. 

Rev. Lett. 39, 737 (1977). 
[115] F. W. Stecker, Phys. Rev. 180, 1264 (1969). 
[116] E. Hayward, Rev. Mod. Phys. 35, 324 (1963). 

[117] M. H. Salamon and F. W. Stecker, Astrophys. J. 493, 547 (1998) [arXiv:astro-ph/9704166]. 

[118] L. N. Epele and E. Roulet, JHEP 9810, 009 (1998) [arXiv:astro-ph/9808104]. 

[119] L. A. Anchordoqui, M. T. Dova, N. Epele and J. D. Swain, in Proceedings of 25th International 

Cosmic Ray Conference (Durban, South Africa) 7, 353 (1997). 
[120] L. A. Anchordoqui, M. T. Dova, L. N. Epele and J. D. Swain, Phys. Rev. D 57, 7103 (1998) 

[arXiv:astro-ph/9708082]. 

[121] L. A. Anchordoqui, J. F. Beacom, H. Goldberg, S. Palomares-Ruiz and T. J. Weiler, Phys. Rev. D 
75, 063001 (2007) [arXiv:astro-ph/0611581]. 



85 



[122] S. Karakula and W. Tkaczyk, Astropart. Phys. 1, 229 (1993). 

[123] F. W. Stecker and M. H. Salamon, Astrophys. J. 512, 521 (1992) [arXiv:astro-ph/9808110]. 
[124] L. A. Anchordoqui, Ph.D. Thesis [arXiv:astro-ph/98 12445]. 

[125] D. Harari, S. Mollerach and E. Roulet, JHEP 08, 022 (1999) [arXiv:astro-ph/9906309]. 
[126] K.-T. Kim, P. P. Kronberg, G. Giovannini and T. Venturi, Nature 341, 720 (1989). 
[127] T. E. Clarke, P. P. Kronberg and H. Bohringer, Astrophys. J. 547, LI 1 1 (2001). 
[128] P. P. Kronberg, Rept. Prog. Phys. 57, 325 (1994). 

[129] J. D. Barrow, P. G. Ferreira and J. Silk, Phys. Rev. Lett. 78, 3610 (1997) [arXiv:astro-ph/9701063]. 
[130] K. Jedamzik, V. Katalinic and A. V. Olinto, Phys. Rev. Lett. 85, 700 (2000) [arXiv:astro- 
ph/9911100]. 

[131] P. Blasi, S. Buries and A. V. Olinto, Astrophys. J. 514, L79 (1999) [arXiv:astro-ph/9812487]. 
[132] G. R. Farrar and T. Piran, Phys. Rev. Lett. 84, 3527 (2000) [arXiv:astro-ph/9906431]. 
[133] G. R. Farrar and T. Piran, arXiv:astro-ph/0010370. 

[134] L. A. Anchordoqui, H. Goldberg and T. J. Weiler, Phys. Rev. Lett. 87, 081 101 (2001) [arXiv:astro- 
ph/0103043]. 

[135] L. A. Anchordoqui, H. Goldberg and T. J. Weiler, arXiv: 1103.0536 [astro-ph.HE]. 
[136] G. Sigl, M. Lemoine and P. Biermann, Astropart. Phys. 10, 141 (1999) [arXiv:astro-ph/9806283]. 
[137] G. F. D. Duff and D. Naylor, Differential Equations of Applied Mathematics, (John Wiley & Sons, 
Inc., New York, 1966). 

[138] M. M. Winn, J. Ulrichs, L. S. Peak, C. B. A. Mccusker and L. Horton, J. Phys. G 12, 653 (1986). 
[139] T. Piran, arXiv: 1005.3311 [astro-ph.HE]. 

[140] M. Lemoine and E. Waxman, JCAP 0911, 009 (2009) [arXiv:0907.1354]. 

[141] M. Ahlers, L. A. Anchordoqui, H. Goldberg, F Halzen, A. Ringwald and T. J. Weiler, Phys. Rev. 

D 72, 023001 (2005) [arXiv:astro-ph/0503229]. 
[142] T. K. Gaisser, arXiv:astro-ph/9707283. 
[143] E. Waxman, Astrophys. J. 452, LI (1995). 

[144] E. Waxman and J. N. Bahcall, Phys. Rev. D 59, 023002 (1999) [arXiv:hep-ph/9807282]. 

[145] J. Alvarez-Muniz, R. Engel, T. K. Gaisser, J. A. Ortiz and T. Stanev, Phys. Rev. D 66, 033011 

(2002) [arXiv:astro-ph/0205302]. 
[146] J. G. Learned and S. Pakvasa, Astropart. Phys. 3, 267 (1995) [arXiv:hep-ph/9405296]. 
[147] E. Waxman and J. N. Bahcall, Phys. Rev. Lett. 78, 2292 (1997) [arXiv:astro-ph/9701231]. 
[148] D. Band et al, Astrophys. J. 413, 281 (1993). 

[149] M. Ahlers, M. C. Gonzalez-Garcia and F Halzen, arXiv: 1103.3421 [astro-ph.HE]. 
[150] L. A. Anchordoqui, D. Hooper, S. Sarkar and A. M. Taylor, Astropart. Phys. 29, 1 (2008) 
[arXiv : astro-ph/070300 1 ] . 

[151] L. A. Anchordoqui, H. Goldberg, D. Hooper, S. Sarkar and A. M. Taylor, Phys. Rev. D 76, 123008 

(2007) [arXiv:0709.0734 [astro-ph]]. 
[152] F. Halzen and A. O'Murchadha, arXiv:0802.0887 [astro-ph]. 

[153] L. A. Anchordoqui, H. Goldberg, F. Halzen and T. J. Weiler, Phys. Lett. B 600, 202 (2004) 

[arXiv:astro-ph/0404387]. 
[154] P. Padovani and C. M. Urry, Astrophys. J. 356, 75 (1990). 
[155] V. S. Beresinsky and G. T. Zatsepin, Phys. Lett. B 28, 423 (1969). 
[156] F. W. Stecker, Astrophys. J. 228, 919 (1979). 

[157] M. Ahlers, L. A. Anchordoqui and S. Sarkar, Phys. Rev. D 79, 083009 (2009) [arXiv:0902.3993]. 
[158] H. Yiiksel, M. D. Kistler, J. F. Beacom and A. M. Hopkins, Astrophys. J. 683, L5 (2008) 



86 



[arXiv:0804.4008 [astro-ph]]. 
[159] A. Franceschini, G. Rodighiero and M. Vaccari, Astron. Astrophys. 487, 837 (2008) 

[arXiv:0805.1841 [astro-ph]]. 
[160] M. Ahlers, L. A. Anchordoqui, M. C. Gonzalez-Garcia, F. Halzen and S. Sarkar, Astropart. Phys. 

34, 106 (2010) [arXiv: 1005.2620 [astro-ph.HE]]. 
[161] F. Aharonian etal. [H.E.S.S. Collaboration], Nature 440, 1018 (2006) [arXiv:astro-ph/0508073]. 
[162] E. Aliu etal. [MAGIC Collaboration], Science 320, 1752 (2008) [arXiv:0807.2822 [astro-ph]]. 
[163] A. A. Abdo et al. [Fermi-LAT Collaboration], Astrophys. J. 723, 1082 (2010) [ arXiv: 1005.0996 

[astro-ph.HE]]. 

[164] A. M. Hopkins and J. F. Beacom, Astrophys. J. 651, 142 (2006) [arXiv:astro-ph/0601463]. 
[165] T. A. Clark, L. W. Brown, and J. K. Alexander, Nature 228, 847 (1970). 

[166] R. J. Protheroe and P. L. Biermann, Astropart. Phys. 6, 45 (1996) [Erratum-ibid. 7, 181 (1996)], 

[arXiv:astro-ph/96051 19]. 
[167] R. J. Protheroe and T. Stanev, Mon. Not. R. Astron. Soc. 264, 191 (1993). 
[168] A. Mucke, R. Engel, J. P. Rachen, R. J. Protheroe and T. Stanev, Comput. Phys. Commun. 124, 

290 (2000) [arXiv:astro-ph/9903478]. 
[169] A. A. Abdo et al. [Fermi-LAT Collaboration], Phys. Rev. Lett. 104, 101101 (2010) 

[arXiv: 1002.3603 [astro-ph.HE]]. 
[170] V. Berezinsky, A. Gazizov, M. Kachelriess and S. Ostapchenko, Phys. Lett. B 695, 13 (2011) 

[arXiv: 1003. 1496 [astro-ph.HE]]. 
[171] H. Yiiksel and M. D. Kistler, Phys. Rev. D 75, 083004 (2007) [arXiv:astro-ph/06 10481]. 
[172] T. Stanev, arXiv:0808.1045 [astro-ph]. 

[173] G. Hasinger, T. Miyaji and M. Schmidt, Astron. Astrophys. 441, 417 (2005) [arXiv:astro- 
ph/0506118]. 

[174] R. Engel, D. Seckel and T. Stanev, Phys. Rev. D 64, 093010 (2001) [arXiv:astro-ph/0101216]. 
[175] L. A. Anchordoqui and T. Montaruli, Ann. Rev. Nucl. Part. Sci. 60, 129 (2010) [arXiv:0912.1035 

[astro-ph.HE]]. 
[176] F. Halzen, Science 315, 66 (2007). 

[177] K. S. Capelle, J. W. Cronin, G. Parente and E. Zas, Astropart. Phys. 8, 321 (1998) [arXiv:astro- 
ph/9801313]. 

[178] X. Bertou, P. Billoir, O. Deligny, C. Lachaud and A. Letessier-Selvon, Astropart. Phys. 17, 183 

(2002) [arXiv:astro-ph/0104452]. 
[179] J. L. Feng, P. Fisher, F. Wilczek and T. M. Yu, Phys. Rev. Lett. 88, 161102 (2002) [arXiv:hep- 

ph/0105067]. 

[180] D. Fargion, Astrophys. J. 570, 909 (2002) [arXiv:astro-ph/0002453]. 

[181] J. Abraham et al. [Pierre Auger Collaboration], Phys. Rev. Lett. 100, 211101 (2008) 

[arXiv:0712.1909 [astro-ph]]. 
[182] J. Abraham et al. [Pierre Auger Collaboration], Phys. Rev. D 79, 102001 (2009) [arXiv:0903.3385 

[astro-ph.HE]]. 

[183] J. Abraham et al. [Pierre Auger Collaboration], arXiv:0906.2347. 
[184] R. Abbasi etal. [IceCube Collaboration], arXiv: 1103.4250 [astro-ph.CO]. 
[185] G. J. Feldman and R. D. Cousins, Phys. Rev. D 57, 3873 (1998) [arXiv:physics/9711021]. 
[186] L. A. Anchordoqui, J. L. Feng, H. Goldberg and A. D. Shapere, Phys. Rev. D 66, 103002 (2002) 
[arXiv:hep-ph/0207139]. 

[187] R. Gandhi, C. Quigg, M. H. Reno and I. Sarcevic, Phys. Rev. D 58, 093009 (1998) [arXiv:hep- 
ph/9807264]. 



87 



[188] A. Capella, U. Sukhatme, C. I. Tan and J. Tran Thanh Van, Phys. Rept. 236, 225 (1994). 
[189] E. Predazzi, arXiv:hep-ph/9809454. 

[190] D. Cline, F. Halzen and J. Luthe, Phys. Rev. Lett. 31, 491 (1973). 

[191] S. D. Ellis and M. B. Kislinger, Phys. Rev. D 9, 2027 (1974). 

[192] F. Halzen, Nucl. Phys. B 92, 404 (1975). 

[193] G. Pancheri and C. Rubbia, Nucl. Phys. A 418, 117C (1984). 

[194] T. K. Gaisser and F. Halzen, Phys. Rev. Lett. 54, 1754 (1985). 

[195] J. Dias de Deus, Nucl. Phys. B 252, 369 (1985). 

[196] G. Pancheri and Y. Srivastava, Phys. Lett. B 159, 69 (1985). 

[197] G. Pancheri and Y. N. Srivastava, Phys. Lett. B 182, 199 (1986). 

[198] C. Albajar etal. [UA1 Collaboration], Nucl. Phys. B 309, 405 (1988). 

[199] C. Adloff etal. [HI Collaboration], Eur. Phys. J. C 21, 33 (2001) [arXiv:hep-ex/00 12053]. 

[200] L. Anchordoqui, M. T. Dova, A. G. Mariazzi, T. McCauley, T C. Paul, S. Reucroft and J. Swain, 

Annals Phys. 314, 145 (2004) [arXiv:hep-ph/0407020]. 
[201] V. N. Gribov and L. N. Lipatov, Yad. Fiz. 15, 1218 (1972) [Sov. J. Nucl. Phys. 15, 675 (1972)]. 
[202] V. N. Gribov and L. N. Lipatov, Yad. Fiz. 15, 781 (1972) [Sov. J. Nucl. Phys. 15, 438 (1972)]. 
[203] Y. L. Dokshitzer, Sov. Phys. JETP 46, 641 (1977) [Zh. Eksp. Teor. Fiz. 73, 1216 (1977)]. 
[204] G. Altarelli and G. Parisi, Nucl. Phys. B 126, 298 (1977). 
[205] M. Dittmar et al, arXiv:0901.2504 [hep-ph]. 

[206] S. Chekanov et al. [ZEUS Collaboration], Phys. Rev. D 67, 012007 (2003) [arXiv:hep- 
ex/0208023]. 

[207] C. Adloff etal. [HI Collaboration], Eur. Phys. J. C 30, 1 (2003) [arXiv:hep-ex/0304003]. 

[208] R. Engel, Nucl. Phys. Proc. Suppl. 122, 40 (2003). 

[209] L. Anchordoqui and F. Halzen, arXiv:0906.1271. 

[210] J. Kwiecinski and A. D. Martin, Phys. Rev. D 43, 1560 (1991). 

[21 1] J. Dias de Deus and J. Kwiecinski, Phys. Lett. B 196, 537 (1987). 

[212] J. Dias De Deus, Nucl. Phys. B 59, 231 (1973). 

[213] U. Amaldi and K. R. Schubert, Nucl. Phys. B 166, 301 (1980). 

[214] R. Castaldi and G. Sanguinetti, Ann. Rev. Nucl. Part. Sci. 35, 351 (1983). 

[215] M. M. Block and R. N. Cahn, Rev. Mod. Phys. 57, 563 (1985). 

[216] R. J. Glauber and G. Matthiae, Nucl. Phys. B 21, 135 (1970). 

[217] P. L'Heureux, B. Margolis and P. Valin, Phys. Rev. D 32, 1681 (1985). 

[218] L. Durand and H. Pi, Phys. Rev. Lett. 58, 303 (1987). 

[219] L. Durand and H. Pi, Phys. Rev. D 38, 78 (1988). 

[220] T. K. Gaisser and T. Stanev, Phys. Lett. B 219, 375 (1989). 

[221] R. S. Fletcher, T. K. Gaisser, P. Lipari and T. Stanev, Phys. Rev. D 50, 5710 (1994). 

[222] N. N. Kalmykov, S. S. Ostapchenko and A. I. Pavlov, Nucl. Phys. Proc. Suppl. 52B, 17 (1997). 

[223] M. Froissart, Phys. Rev. 123, 1053 (1961). 

[224] J. Ranft, Phys. Rev. D 51, 64 (1995). 

[225] E. J. Ahn, R. Engel, T K. Gaisser, P. Lipari and T. Stanev, Phys. Rev. D 80, 094003 (2009) 

[arXiv:0906.4113 [hep-ph]]. 
[226] S. Ostapchenko, AIP Conf. Proc. 928, 118 (2007) [arXiv:0706.3784 [hep-ph]]. 
[227] S. Roesler, R. Engel and J. Ranft, arXiv:hep-ph/0012252. 
[228] T. Sjostrand, Int. J. Mod. Phys. A 3, 751 (1988). 



88 



[229] J. Engel, T. K. Gaisser, T. Stanev and P. Lipari, Phys. Rev. D 46, 5013 (1992). 
[230] R. D. Woods and D. S. Saxon, Phys. Rev. 95, 577 (1954). 

[231] L. Portugal and T. Kodama, Nucl. Phys. A 837, 1 (2010) [arXiv:0910.3701 [hep-ph]]. 
[232] K. Belov [HiRes Collaboration], Nucl. Phys. Proc. Suppl. 151, 197 (2006). 
[233] M. M. Block, F. Halzen and T. Stanev, Phys. Rev. Lett. 83, 4926 (1999) [arXiv:hep-ph/9908222]. 
[234] M. M. Block, Phys. Rev. D 76, 111503 (2007) [arXiv:0705.3037 [hep-ph]]. 
[235] G. Aielli et al [ARGO-YBJ Collaboration], Phys. Rev. D 80, 092004 (2009) [arXiv:0904.4198 
[hep-ex]]. 

[236] M. Aglietta et al, Phys. Rev. D 79, 032004 (2009). 

[237] H. H. Mielke, M. Foeller, J. Engler and J. Knapp, J. Phys. G 20, 637 (1994). 
[238] M. Honda et al, Phys. Rev. Lett. 70, 525 (1993). 

[239] M. N. Dyakonov et al, in Proceedings of 21st International Cosmic Ray Conference (Adelaide, 
Australia) 9, 252 (1990). 

[240] S. P. Knurenko, V. R. Sleptsova, I. E. Sleptsov, N. N. Kalmykov and S. S. Ostapchenko, in Pro- 
ceedings of 26th International Cosmic Ray Conference (Salt Lake City, Utah) 1, 372 (1999). 

[241] R. A. Nam, S. I. Nikolsky, V. P. Pavlyuchenko, A. P. Chubenko and V. I. Yakovlev, in Proceedings 
of 14th International Cosmic Ray Conference (Munich, Germany) 7, 2258 (1975). 

[242] R. M. Baltrusaitis et al, Phys. Rev. Lett. 52, 1380 (1984). 

[243] T. Antoni etal [KASCADE Collaboration], I. Phys. G 27, 1785 (2001) [arXiv:astro-ph/0 106494]. 
[244] K. Eggert, Nucl. Phys. Proc. Suppl. 122, 447 (2003). 

[245] ATLAS Collaboration, Measurement of the inelastic proton-proton cross section at = 7 TeV 

with the ATLAS detector, ATLAS-CONF-201 1-002 (201 1). 
[246] A. Achilli, R. M. Godbole, A. Grau, G Pancheri, O. Shekhovtsova and Y. N. Srivastava, 

arXiv: 1102. 1949 [hep-ph]. 
[247] M. M. Block and F. Halzen, arXiv: 1 102.3 163 [hep-ph]. 
[248] J. Alcaraz, these proceedings. 

[249] M. Gluck, E. Reya and A. Vogt, Z. Phys. C 53, 127 (1992). 
[250] M. Gluck, E. Reya and A. Vogt, Z. Phys. C 67, 433 (1995). 

[251] M. Gluck, E. Reya and A. Vogt, Eur. Phys. I. C 5, 461 (1998) [arXiv:hep-ph/9806404]. 
[252] A. D. Martin, R. G. Roberts, W. J. Stirling and R. S. Thorne, Eur. Phys. J. C 4, 463 (1998) 
[arXiv:hep-ph/9803445]. 

[253] H. L. Lai et al [CTEQ Collaboration], Eur. Phys. J. C 12, 375 (2000) [arXiv:hep-ph/9903282]. 
[254] A. Achilli, R. Hegde, R. M. Godbole, A. Grau, G. Pancheri and Y. Srivastava, Phys. Lett. B 659, 

137 (2008) [arXiv:0708.3626 [hep-ph]]. 
[255] F. Bloch and A. Nordsieck, Phys. Rev. 52, 54 (1937). 

[256] A. Grau, R. M. Godbole, G Pancheri and Y. N. Srivastava, Phys. Lett. B 682, 55 (2009) 

[arXiv:0908.1426 [hep-ph]]. 
[257] A. Grau, G. Pancheri, O. Shekhovtsova and Y. N. Srivastava, Phys. Lett. B 693, 456 (2010) 

[arXiv: 1008.41 19 [hep-ph]]. 
[258] M. M. Block and F. Halzen, Phys. Rev. D 72, 036006 (2005) [Erratum-ibid. D 72, 039902 (2005)] 

[arXiv:hep-ph/0506031]. 
[259] M. M. Block and F. Halzen, Phys. Rev. D 73, 054022 (2006) [arXiv:hep-ph/0510238]. 
[260] M. M. Block, E. M. Gregores, F. Halzen and G. Pancheri, Phys. Rev. D 60, 054024 (1999) 

[arXiv:hep-ph/9809403]. 
[261] M. M. Block, Phys. Rept. 436, 71 (2006) [arXiv:hep-ph/0606215]. 



89 



[262] D. d'Enterria, R. Engel, T. Pierog, S. Ostapchenko and K. Werner, arXiv: 1101.5596 [astro-ph.HE]. 
[263] H. M. J. Barbosa, F. Catalani, J. A. Chinellato and C. Dobrigkeit, arXiv:astro-ph/03 10234. 
[264] T. J. L. Mccomb, R. J. Protheroe and K. E. Turver, J. Phys. G 5, 1613 (1979). 
[265] B. Rossi and K. Greisen, Rev. Mod. Phys. 13, 240 (1941). 
[266] H. Bethe and W. Heitler, Proc. Roy. Soc. Lond. A 146, 83 (1934). 
[267] S. Eidelman et al. [Particle Data Group], Phys. Lett. B 592, 1 (2004). 
[268] B. Rossi, High Energy Particles (Prentice-Hall, Inc., Englewood Cliffs, NY, 1952). 
[269] R. C. Weast, CRC Handbook of Chemistry and Physics, (CRC Press, Boca Raton, FL, USA, 
1981). 

[270] Y. S. Tsai, Rev. Mod. Phys. 46, 815 (1974) [Erratum Rev. Mod. Phys. 49, 421 (1977)]. 
[271] T. Stanev, C. Vankov, R. E. Streitmatter, R. W. Ellsworth and T. Bowen, Phys. Rev. D 25, 1291 
(1982). 

[272] A. N. Cillis, H. Fanchiotti, C. A. Garcia Canal and S. J. Sciutto, Phys. Rev. D 59, 113012 (1999) 

[arXiv:astro-ph/9809334]. 
[273] S. J. Sciutto, arXiv:astro-ph/99 11331. 

[274] H. P. Vankov, N. Inoue and K. Shinozaki, Phys. Rev. D 67, 043002 (2003) [arXiv:astro- 
ph/0211051]. 

[275] X. Bertou, P. Billoir and S. Dagoret-Campagne, Astropart. Phys. 14, 121 (2000). 
[276] L. Anchordoqui and H. Goldberg, Phys. Lett. B 583, 213 (2004) [arXiv:hep-ph/03 10054]. 
[277] A. N. Cillis and S. J. Sciutto, Phys. Rev. D 64, 013010 (2001) [arXiv:astro-ph/0010488]. 
[278] W. Heitler. The Quantum Theory of Radiation, 2nd. Edition, (Oxford University Press, London, 
1944). 

[279] J. Linsley, in Proceedings of 15th International Cosmic Ray Conference (Plovdiv, Bulgaria) 12, 
89 (1977). 

[280] J. Matthews, Astropart. Phys. 22, 387 (2005). 

[281] R. Ulrich, J. Blumer, R. Engel, F. Schussler and M. Unger, New J. Phys. 11, 065018 (2009) 

[arXiv:0903.0404 [astro-ph.HE]]. 
[282] M. Ave, R. A. Vazquez, E. Zas, J. A. Hinton and A. A. Watson, Astropart. Phys. 14, 109 (2000) 

[arXiv:astro-ph/000301 1]. 
[283] A. V. Plyasheshnikov and F. A. Aharonian, J. Phys. G 28, 267 (2002) [arXiv:astro-ph/0107592]. 
[284] D. Heck, G. Schatz, T. Thouw, J. Knapp and J. N. Capdevielle, FZKA-6019 (1998). 
[285] J. Knapp, D. Heck, S. J. Sciutto, M. T. Dova and M. Risse, Astropart. Phys. 19, 77 (2003) 

[arXiv: astro-ph/02064 1 4] . 
[286] A. Kusenko and T. J. Weiler, Phys. Rev. Lett. 88, 161101 (2002) [arXiv:hep-ph/0 106071]. 
[287] L. A. Anchordoqui, J. L. Feng, H. Goldberg and A. D. Shapere, Phys. Rev. D 65, 124027 (2002) 

[arXiv:hep-ph/01 12247]. 
[288] D. Gora, M. Roth and A. Tamburro, Astropart. Phys. 26, 402 (2007). 

[289] A. Gazizov and M. P. Kowalski, Comput. Phys. Commun. 172, 203 (2005) [arXiv:astro- 
ph/0406439]. 

[290] A. Cooper-Sarkar and S. Sarkar, JHEP 0801, 075 (2008) [arXiv:07 10.5303 [hep-ph]]. 
[291] T. Sjostrand, S. Mrenna and P. Z. Skands, JHEP 0605, 026 (2006) [arXiv:hep-ph/0603175]. 
[292] S. Jadach, Z. Was, R. Decker and J. H. Kuhn, Comput. Phys. Commun. 76, 361 (1993). 
[293] Consortium for Spatial Information (CGIAR-CSI). http : //srtm . csi . cgiar . org/ 
[294] L. A. Anchordoqui, H. Goldberg, D. Gora, T. Paul, M. Roth, S. Sarkar and L. L. Winders, Phys. 
Rev. D 82, 043001 (2010) [arXiv: 1004.3 190 [hep-ph]]. 



90 



[295] S. J. Sciutto, arXiv:astro-ph/0106044. 

[296] S. Argiro et al, Nucl. Instrum. Meth. A 580, 1485 (2007) [arXiv:0707.1652 [astro-ph]]. 
[297] M. Gluck, S. Kretzer and E. Reya, Astropart. Phys. 11, 327 (1999) [arXiv:astro-ph/9809273]. 
[298] P. M. Nadolsky et al, Phys. Rev. D 78, 013004 (2008) [arXiv:0802.0007 [hep-ph]]. 
[299] N. Armesto, C. Merino, G. Parente and E. Zas, Phys. Rev. D 77, 013001 (2008) [arXiv: 0709. 4461 
[hep-ph]]. 

[300] E. V. Bugaev and Yu. V. Shlepin, Phys. Rev. D 67, 034027 (2003) [arXiv:hep-ph/0203096]. 
[301] H. Abramowicz and A. Levy, arXiv:hep-ph/9712415. 

[302] A. Capella, A. Kaidalov, C. Merino and J. Tran Thanh Van, Phys. Lett. B 337, 358 (1994) 
[arXiv:hep-ph/9405338]. 

[303] S. I. Dutta, M. H. Reno, I. Sarcevic and D. Seckel, Phys. Rev. D 63, 094020 (2001) [arXiv:hep- 

ph/00 12350]. 
[304] G. Cowan, these proceedings. 

[305] S. Baker and R. D. Cousins, Nucl. Instrum. Meth. A 221, 437 (1984). 

[306] L. A. Anchordoqui, A. M. Cooper-Sarkar, D. Hooper and S. Sarkar, Phys. Rev. D 74, 043008 

(2006) [arXiv:hep-ph/0605086]. 
[307] T. J. Weiler, Nucl. Phys. Proc. Suppl. 134, 47 (2004). 

[308] L. A. Anchordoqui, Z. Fodor, S. D. Katz, A. Ringwald and H. Tu, JCAP 0506, 013 (2005) 
[arXiv:hep-ph/0410136]. 

[309] L. A. Anchordoqui, H. Goldberg, F. Halzen and T. J. Weiler, Phys. Lett. B 593, 42 (2004) 

[arXiv: astro-ph/03 1 1 002] . 
[310] D. Hooper, A. Taylor and S. Sarkar, Astropart. Phys. 23, 11 (2005) [arXiv:astro-ph/0407618]. 
[311] M. Ave, N. Busca, A. V. Olinto, A. A. Watson and T. Yamamoto, Astropart. Phys. 23, 19 (2005) 

[arXiv: astro-ph/04093 1 6] . 

[312] L. A. Anchordoqui, J. F. Beacom, H. Goldberg, S. Palomares-Ruiz and T. J. Weiler, Phys. Rev. 

Lett. 98, 121101 (2007) [arXiv:astro-ph/06 11580]. 
[313] X. Y. Wang, R. Y. Liu and F. Aharonian, arXiv: 1103.3574 [astro-ph.HE]. 

[314] L. A. Anchordoqui, H. Goldberg, F. Halzen and T. J. Weiler, Phys. Lett. B 621, 18 (2005) 

[arXiv:hep-ph/0410003]. 
[315] G. R. Blumenthal and R. J. Gould, Rev. Mod. Phys. 42, 237 (1970). 
[316] S. Lee, Phys. Rev. D 58, 043004 (1998) [arXiv:astro-ph/9604098]. 

[317] S. V. Demidov and O. E. Kalashev, J. Exp. Theor. Phys. 108, 764 (2009) [arXiv:08 12.0859 [astro- 
ph]]. 

[318] J. Wdowczyk, W. Tkaczyk and A. W. Wolfendale, J. Phys. A 5, 1419 (1972). 

[319] A. Neronov and I. Vovk, Science 328, 73 (2010) [arXiv: 1006.3504 [astro-ph.HE]]. 

[320] A. M. Taylor, I. Vovk and A. Neronov, arXiv: 1101. 0932 [astro-ph.HE]. 

[321] G. Matthiae, Nucl. Phys. Proc. Suppl. 99A, 281 (2001). 

[322] P. Lipari, Nucl. Phys. Proc. Suppl. 122, 133 (2003) [arXiv:astro-ph/0301196]. 

[323] C. Augier et al [UA4/2 Collaboration], Phys. Lett. B 315, 503 (1993). 

[324] TOTEM Technical Design Report, CERN-LHCC-2004-002. 

[325] G. 't Hooft, Phys. Rev. Lett. 37, 8 (1976). 

[326] G. 't Hooft, Phys. Rev. D 14, 3432 (1976) [Erratum-ibid. D 18, 2199 (1978)]. 
[327] F. R. Klinkhamer and N. S. Manton, Phys. Rev. D 30, 2212 (1984). 
[328] V. A. Kuzmin, V. A. Rubakov and M. E. Shaposhnikov, Phys. Lett. B 155, 36 (1985). 
[329] M. Fukugita and T. Yanagida, Phys. Lett. B 174, 45 (1986). 



91 



[330] P. B. Arnold and L. D. McLerran, Phys. Rev. D 36, 581 (1987). 

[331] H. Aoyama and H. Goldberg, Phys. Lett. B 188, 506 (1987). 

[332] A. Ringwald, Nucl. Phys. B 330, 1 (1990). 

[333] O. Espinosa, Nucl. Phys. B 343, 310 (1990). 

[334] A. Ringwald, JHEP 0310, 008 (2003) [arXiv:hep-ph/0307034]. 

[335] F. L. Bezrukov, D. Levkov, C. Rebbi, V. A. Rubakov and P. Tinyakov, Phys. Rev. D 68, 036005 

(2003) [arXiv:hep-ph/0304180]. 
[336] F. L. Bezrukov, D. Levkov, C. Rebbi, V. A. Rubakov and P. Tinyakov, Phys. Lett. B 574, 75 (2003) 

[arXiv:hep-ph/0305300]. 
[337] T. M. Gould and S. D. H. Hsu, Mod. Phys. Lett. A 9, 1589 (1994) [arXiv:hep-ph/931 1291]. 



92 



